NWGOM Population Genetics
Ryan Eckert
version: 04 October, 2024
About this document
All analyses preformed with R version 4.4.1.
Basic setup of R environment
Loading required packages
For the following analyses we will require the use of a number of different R packages. We can use the following code to quickly load in the packages and install any packages not previously installed in the R console.
if (!require("pacman")) install.packages("pacman")
pacman::p_load_gh("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis", "ropensci/rnaturalearthhires", "KarstensLab/microshades")
pacman::p_load("cowplot", "car", "ggrepel", "ggspatial", "paletteer", "patchwork", "rnaturalearth", "sf", "Hmisc", "MCMC.OTU", "pairwiseAdonis", "RColorBrewer", "Redmonder", "flextable", "lubridate", "officer", "adegenet", "dendextend", "gdata", "ggdendro", "hierfstat", "kableExtra", "poppr", "reshape2", "StAMPP", "vcfR", "vegan", "boa", "magick", "sdmpredictors", "ggcorrplot", "tidyverse", "TeachingDemos", "LaplacesDemon", "adespatial", "ggnewscale", "ggbeeswarm", "multcomp", "rstatix", "R.utils", "graph4lg", "ggstance", "ggtreeExtra", "ggtree")
options("scipen" = 10)
# load("nwgom.RData")Themes and colors
Making color palettes to use throughout all plots
gomPal = c("#D53E4F", "#F88D51", "#FEE08B", "#66C2A5", "#3288BD")
boundPal = c("gray30", paletteer_d("vapoRwave::vapoRwave")[10])
pink = "#FF6A8BFF"
purple = paletteer_d("vapoRwave::vapoRwave")[10]
# mcavKColPal = c(colorRampPalette(colors = c("#1B9E77","#FFFFFF"))(3)[1:2], "azure3")
# sintKColPal = c(colorRampPalette(colors = c("#D95F02","#FFFFFF"))(3)[1:2], "azure3")
# ofavKColPal = colorRampPalette(colors = c("#7570B3","#FFFFFF"))(4)[1:3]
# xestoKColPal = c(colorRampPalette(colors = c("#E7298A","#FFFFFF"))(7)[c(1:6)], "azure3")
mcavKColPal = c(colorRampPalette(colors = c("#1B9E77","#FFFFFF"))(7)[c(1,5,2,3,4,6)], "azure3")
sintKColPal = c(colorRampPalette(colors = c("#D95F02","#FFFFFF"))(7)[c(1,5,2,3,4,6)], "azure3")
ofavKColPal = colorRampPalette(colors = c("#7570B3","#FFFFFF"))(7)[c(1,5,3,2,4,6)]
xestoKColPal = c(colorRampPalette(colors = c("#E7298A","#FFFFFF"))(7)[c(1:6)], "azure3")
profPal = rev(c(microshades_palette("micro_green", 5), microshades_palette("micro_cvd_turquoise", 5), microshades_palette("micro_cvd_orange", 3),microshades_palette("micro_cvd_purple", 1, lightest = F), microshades_palette("micro_purple", 5)))
colPalZoox = c("#6A4C93", "#1BE7FF", "#C6FCA2", "#FF6A8B")Plot themes
dendTheme = theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 10, color = "black", angle = 90),
axis.text.y = element_text(size = 8, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
legend.key.size = unit(0.75, "line"),
legend.key = element_blank(),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
legend.direction = "horizontal",
legend.box = "vertical",
legend.position = c(0.13, 0.1))
pcaTheme = theme(axis.title.x = element_text(color = "black", size = 10),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.line.x = element_blank(),
axis.title.y = element_text(color = "black", size = 10),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.y = element_blank(),
legend.position = "none",
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
legend.key.size = unit(5, "pt"),
panel.border = element_rect(color = "black", linewidth = 1),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
admixTheme = theme(
panel.grid = element_blank(),
panel.background = element_rect(fill = "gray85"),
plot.background = element_blank(),
panel.border = element_rect(fill = NA, color = "black", linewidth = 0.75, linetype = "solid"),
panel.spacing.x = grid:::unit(0.05, "lines"),
panel.spacing.y = grid:::unit(0.05, "lines"),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.title = element_blank(),
strip.background.x = element_blank(),
strip.background.y = element_blank(),
strip.text = element_text(size = 8),
strip.text.y.left = element_text(size = 10, angle = 90),
strip.text.x.bottom = element_text(vjust = 1, color = NA),
legend.key = element_blank(),
legend.position = "none",
legend.title = element_text(size = 8))
violinTheme = theme(axis.title = element_text(color = "black", size = 12),
axis.ticks.x = element_blank(),
axis.text.x = element_text(color = "black", size = 10),
legend.position = "none",
legend.key.size = unit(0.3, 'cm'),
panel.border = element_rect(color = "black", size = 1),
panel.background = element_blank(),
plot.background = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
zooxTheme = theme(plot.title = element_text(),
panel.grid = element_blank(),
# panel.background = element_blank(),
panel.background = element_rect(fill = "gray85"),
panel.border = element_rect(fill = NA, color = "black", size = 0.75, linetype = "solid"),
panel.spacing.x = grid:::unit(0.05, "lines"),
panel.spacing.y = grid:::unit(0.82, "lines"),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks.x = element_blank(),
axis.ticks.y = element_blank(),
axis.title = element_blank(),
strip.background.x = element_blank(),
strip.background.y = element_blank(),
strip.text = element_text(size = 12),
strip.text.y.left = element_text(size = 12, angle = 90),
strip.text.x.bottom = element_text(vjust = -.1, color = NA),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
legend.key.size = unit(0.75, "line"),
legend.key = element_blank(),
legend.position = "right",
legend.direction = "vertical",
legend.box = "vertical")Sampling info
Map of study sites
nwgomSites = read.csv("../data/nwgomMetaData.csv", header = TRUE)
nwgomSites$depthZone = factor(nwgomSites$depthZone)
nwgomSites$depthZone = factor(nwgomSites$depthZone, levels = levels(nwgomSites$depthZone)[c(2,1)])
nwgomSites$site = factor(nwgomSites$site)
nwgomSites$bank = factor(nwgomSites$bank)
nwgomSites$bank = factor(nwgomSites$bank, levels = levels(nwgomSites$bank)[c(5, 2, 1, 3, 4)])
nwgomSites$date = mdy(nwgomSites$date) %>% format("%d %b %Y")
nwgomBanks = nwgomSites %>% group_by(bank) %>% summarise(latDD = mean(latDD), longDD = mean(longDD), n = n()) %>% droplevels()
nwgomSampleSites = nwgomSites %>% group_by(bank, site, depthZone) %>% summarise(latDD = min(latDD), longDD = min(longDD))## `summarise()` has grouped output by 'bank', 'site'. You can override using the
## `.groups` argument.fgbnmsBounds = read_sf("../data/shp/FGBNMS_py.shp") %>% st_transform(crs = 4326)
fgbnmsBounds$Bank = c("NA","NA","Geyer","NA","Bright","WFGB","NA","McGrail","NA","NA","NA","McGrail","EFGB","NA","NA","NA","NA","NA","NA")
fgbnmsBounds$Bank = factor(fgbnmsBounds$Bank)
fgbnmsBounds$Bank = factor(fgbnmsBounds$Bank, levels = levels(fgbnmsBounds$Bank)[c(6, 2, 1, 3, 4,5)])
states = st_as_sf(ne_states(country = c("United States of America")), scale = "count", crs = 4326) %>% filter(name_en %in% c("Florida", "Alabama","Mississippi", "Louisiana", "Texas", "Georgia", "South Carolina"))
countries = st_as_sf(ne_countries(country = c("Cuba", "Mexico", "Cayman Islands", "The Bahamas", "Belize", "Guatemala"), scale = "Large"), crs = 4326)
# bahamas = read_sf("../data/shp/bahamasShoreline.shp") %>% st_transform(crs = 4326)
# cuba = read_sf("../data/shp/cubaShoreline.shp") %>% st_transform(crs = 4326)
# florida = read_sf("../data/shp/floridaShoreline.shp") %>% st_transform(crs = 4326)
wfgbBathy = read_sf("../data/shp/wb_10m_Cleaned.shp") %>% st_transform(crs = 4326)
efgbBathy = read_sf("../data/shp/eb_10m_Cleaned.shp") %>% st_transform(crs = 4326)
brightBathy = read_sf("../data/shp/Bright_10m_Cleaned.shp") %>% st_transform(crs = 4326)
geyerBathy = read_sf("../data/shp/Geyer_10m_Cleaned.shp") %>% st_transform(crs = 4326)
mcgrailBathy = read_sf("../data/shp/McGrail_10m_Cleaned.shp") %>% st_transform(crs = 4326)nwgomMap = ggplot() +
geom_sf(data = fgbnmsBounds %>% filter(Bank != "NA"), aes(fill = Bank), linewidth = 0.25, color = "black") +
geom_sf(data = fgbnmsBounds %>% filter(Bank == "NA"), fill = "black", linewidth = 0.25, color = "black", alpha = 0.1) +
geom_sf(data = states, fill = "white", linewidth = 0.3) +
scale_fill_manual(values = c(gomPal,"black"), name = "Bank",) +
geom_point(data = nwgomSites, aes(x = longDD, y = latDD, shape = depthZone), color = NA, fill = NA) +
scale_shape_manual(values = c(21, 23), name = "Depth") +
coord_sf(xlim = c(-94.5, -91.5), ylim = c(26.75, 29.75)) +
annotation_scale(location = "br") +
guides(fill = guide_legend(override.aes = list(shape = 22, color = "black", size = 2.5, stroke = 0.25), order = 1), shape = guide_legend(override.aes = list(size = c(2.5, 2.25), stroke = 0.25, color = "black"), order = 2), color = "none") +
theme_bw() +
theme(panel.background = element_rect(fill = "aliceblue"),
plot.background = element_blank(),
panel.border = element_rect(color = "black", size = 1, fill = NA),
axis.title = element_blank(),
axis.ticks = element_line(color = "black"),
axis.text = element_text(color = "black"),
legend.position = c(0.9085, 0.882),
legend.box.background = element_rect(linewidth = 0.45, fill = "white"),
legend.title = element_text(color = "black", size = 8),
legend.text = element_text(color = "black", size = 8),
legend.spacing = unit(-5, "pt"),
legend.key.size = unit(5, "pt"),
legend.key = element_blank(),
legend.background = element_blank()
)
nwgomMap
largeMap = ggplot() +
geom_sf(data = states, fill = "white", linewidth = 0.3) +
geom_sf(data = countries, fill = "white", linewidth = 0.3) +
geom_sf(data = fgbnmsBounds, fill = "black", linewidth = 0.1, color = "black", alpha = 0.1) +
geom_rect(aes(xmin = -95, xmax = -91, ymin = 26, ymax = 30), fill = NA, color = "black", linewidth = 0.75) +
coord_sf(xlim = c(-99, -80), ylim = c(18, 32)) +
theme_bw() +
theme(legend.title = element_text(size = 9, face = "bold"),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
panel.background = element_rect(fill = "aliceblue"),
panel.border = element_rect(color = "black", size = 1, fill = NA),
legend.position = "none",
plot.background = element_blank())
largeMap
inset = ggplot() +
geom_sf(data = wfgbBathy, color = "gray75", size = 0.25) +
geom_sf(data = efgbBathy, color = "gray75", size = 0.25) +
geom_sf(data = brightBathy, color = "gray75", size = 0.25) +
geom_sf(data = geyerBathy, color = "gray75", size = 0.25) +
geom_sf(data = mcgrailBathy, color = "gray75", size = 0.25) +
geom_sf(data = fgbnmsBounds %>% filter(Bank != "NA"), aes(fill = Bank), linewidth = 0.4, color = "black", alpha = 0.2) +
geom_sf(data = fgbnmsBounds %>% filter(Bank == "NA"), fill = "black", linewidth = 0.4, color = "black", alpha = 0.2) +
geom_point(data = nwgomSampleSites, aes(x = longDD, y = latDD, fill = bank, shape = depthZone, size = depthZone), alpha = 0.75) +
scale_fill_manual(values = gomPal, name = "Bank") +
scale_shape_manual(values = c(21, 23), name = "Depth") +
scale_size_manual(values = c(2.5, 2.25)) +
theme_bw() +
theme(legend.title = element_text(size = 9, face = "bold"),
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
panel.background = element_rect(fill = "aliceblue"),
panel.border = element_rect(color = "black", size = 1, fill = NA),
legend.position = "none",
plot.background = element_blank())
wfgb = inset +
geom_label(aes(x = -93.8269, y = 27.882), label = "WFGB", fill = gomPal[1], color = "black", size = 3) +
annotation_scale(location = "br") +
coord_sf(xlim = c(-93.83, -93.808), ylim = c(27.867, 27.883)) +
scale_x_continuous(n.breaks = 10) +
scale_y_continuous(n.breaks = 10)
efgb = inset +
geom_label(aes(x = -93.6074, y = 27.92435), label = "EFGB", fill = gomPal[2], color = "black", size = 3) +
annotation_scale(location = "br") +
coord_sf(xlim = c(-93.61, -93.59), ylim = c(27.905, 27.925)) +
scale_x_continuous(n.breaks = 10) +
scale_y_continuous(n.breaks = 10)
bright = inset +
geom_label(aes(x = -93.3135, y = 27.8984), label = "Bright", fill = gomPal[3], color = "black", size = 3) +
annotation_scale(location = "br") +
coord_sf(xlim = c(-93.316, -93.296), ylim = c(27.879, 27.899)) +
scale_x_continuous(n.breaks = 10) +
scale_y_continuous(n.breaks = 10)
geyer = inset +
geom_label(aes(x = -93.07329, y = 27.8614), label = "Geyer", fill = gomPal[4], color = "black", size = 3) +
annotation_scale(location = "br") +
coord_sf(xlim = c(-93.076, -93.056), ylim = c(27.842, 27.862)) +
scale_x_continuous(n.breaks = 10) +
scale_y_continuous(n.breaks = 10)
mcgrail = inset +
geom_label(aes(x = -92.6018, y = 27.9691), label = "McGrail", fill = gomPal[5], color = "black", size = 3) +
annotation_scale(location = "br") +
coord_sf(xlim = c(-92.605, -92.585), ylim = c(27.955, 27.97)) +
scale_x_continuous(n.breaks = 10) +
scale_y_continuous(n.breaks = 10)
map = (nwgomMap + inset_element(largeMap, top = 1.12, right = 0.33, bottom = 0.63, left = -0.0075, ignore_tag = TRUE) +
inset_element(wfgb, top = 0.8, right = 0.2875, bottom = 0.5, left = -0.0075, ignore_tag = TRUE) +
inset_element(efgb, top = 0.33, right = 0.2875, bottom = 0.04, left = -0.0075, ignore_tag = TRUE) +
inset_element(bright, top = 0.33, right = 0.68, bottom = 0.04, left = 0.33, ignore_tag = TRUE) +
inset_element(geyer, top = 0.33, right = 1.0, bottom = 0.04, left = 0.705, ignore_tag = TRUE) +
inset_element(mcgrail, top = 0.8, right = 1.0, bottom = 0.5, left = 0.7, ignore_tag = TRUE)
)
ggsave("../figures/figure1.png", plot = map, height = 7, width = 7, units = "in", dpi = 300)Detecting clones
Montastraea cavernosa
mcavCloneBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "M. cavernosa") %>% dplyr::filter(!Sample %in% c("MGM010")) # list of bam files
mcavCloneMa = as.matrix(read.table("../data/mcav/mcavClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(mcavCloneMa) = list(mcavCloneBams[,1],mcavCloneBams[,1])
mcavClonesHc = hclust(as.dist(mcavCloneMa),"ave")
mcavClonePops = mcavCloneBams$bank
mcavCloneDepth = mcavCloneBams$depthZone
mcavCloneDend = mcavCloneMa %>% as.dist() %>% hclust(.,"ave") %>% as.dendrogram()
mcavCloneDData = mcavCloneDend %>% dendro_data()
# Making the branches hang shorter so we can easily see clonal groups
mcavCloneDData$segments$yend2 = mcavCloneDData$segments$yend
for(i in 1:nrow(mcavCloneDData$segments)) {
if (mcavCloneDData$segments$yend2[i] == 0) {
mcavCloneDData$segments$yend2[i] = (mcavCloneDData$segments$y[i] - 0.01)}}
mcavCloneDendPoints = mcavCloneDData$labels
mcavCloneDendPoints$pop = mcavClonePops[order.dendrogram(mcavCloneDend)]
mcavCloneDendPoints$depth=mcavCloneDepth[order.dendrogram(mcavCloneDend)]
rownames(mcavCloneDendPoints) = mcavCloneDendPoints$label
# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(mcavCloneDData$segments)) {
if (mcavCloneDData$segments$yend[i] == 0) {
point[i] = mcavCloneDData$segments$y[i] - 0.01
} else {
point[i] = NA}}
mcavCloneDendPoints$y = point[!is.na(point)]
mcavTechReps = c("MGM002.1", "MGM002.2", "MGM002.3", "MGM008.1", "MGM008.2", "MGM008.3", "MGM013.1", "MGM013.2", "MGM013.3", "MGM024.1", "MGM024.2", "MGM038.1", "MGM038.2", "MGM038.3", "MGM072.1", "MGM072.2", "MGM072.3")
mcavCloneDendPoints$depth = factor(mcavCloneDendPoints$depth)
mcavCloneDendPoints$depth = factor(mcavCloneDendPoints$depth,levels(mcavCloneDendPoints$depth)[c(2,1)])
mcavCloneDendPoints$pop = factor(mcavCloneDendPoints$pop)
mcavCloneDendPoints$pop = factor(mcavCloneDendPoints$pop,levels(mcavCloneDendPoints$pop)[c(5, 2, 1, 4, 3)])
mcavCloneDendA = ggplot() +
geom_segment(data = segment(mcavCloneDData), aes(x = x, y = y, xend = xend, yend = yend2), size = 0.5) +
geom_point(data = mcavCloneDendPoints, aes(x = x, y = y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
scale_fill_manual(values =gomPal, name= "Bank")+
scale_shape_manual(values = c(21, 23), name = "Depth Zone")+
geom_hline(yintercept = 0.1575, color = pink, lty = 5, size = 0.75) + # creating a dashed line to indicate a clonal distance threshold
geom_text(data = subset(mcavCloneDendPoints, subset = label %in% mcavTechReps), aes(x = x, y = (y - .015), label = label), angle = 90) + # spacing technical replicates further from leaf
geom_text(data = subset(mcavCloneDendPoints, subset = !label %in% mcavTechReps), aes(x = x, y = (y - .010), label = label), angle = 90) +
labs(y = "Genetic distance (1 - IBS)") +
guides(fill = guide_legend(override.aes = list(shape = 22)))+
coord_cartesian(xlim = c(0,182), ylim = c(0.075, 0.3),expand = FALSE) +
theme_classic()
mcavCloneDend = mcavCloneDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 12, color = "black", angle = 90),
axis.text.y = element_text(size = 10, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "bottom")
mcavCloneDend
ggsave("../figures/extras/mcavCloneDend.png", plot = mcavCloneDend, height = 9, width = 32, units = "in", dpi = 300)Stephanocoenia intersepta
sintCloneBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "S. intersepta") %>% dplyr::filter(!Sample %in% c("SGM013", "SGM169", "SGM177"))
sintCloneMa = as.matrix(read.table("../data/sint/sintClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(sintCloneMa) = list(sintCloneBams[,1],sintCloneBams[,1])
sintClonesHc = hclust(as.dist(sintCloneMa),"ave")
sintClonePops = sintCloneBams$bank
sintCloneDepth = sintCloneBams$depthZone
sintCloneDend = sintCloneMa %>% as.dist() %>% hclust(.,"ave") %>% as.dendrogram()
sintCloneDData = sintCloneDend %>% dendro_data()
# Making the branches hang shorter so we can easily see clonal groups
sintCloneDData$segments$yend2 = sintCloneDData$segments$yend
for(i in 1:nrow(sintCloneDData$segments)) {
if (sintCloneDData$segments$yend2[i] == 0) {
sintCloneDData$segments$yend2[i] = (sintCloneDData$segments$y[i] - 0.01)}}
sintCloneDendPoints = sintCloneDData$labels
sintCloneDendPoints$pop = sintClonePops[order.dendrogram(sintCloneDend)]
sintCloneDendPoints$depth=sintCloneDepth[order.dendrogram(sintCloneDend)]
rownames(sintCloneDendPoints) = sintCloneDendPoints$label
# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(sintCloneDData$segments)) {
if (sintCloneDData$segments$yend[i] == 0) {
point[i] = sintCloneDData$segments$y[i] - 0.01
} else {
point[i] = NA}}
sintCloneDendPoints$y = point[!is.na(point)]
sintTechReps = c("SGM027.1", "SGM027.2", "SGM027.3", "SGM046.1", "SGM046.2", "SGM046.3", "SGM152.1", "SGM152.2", "SGM152.3", "SGM186.1", "SGM186.2", "SGM186.3", "SGM197.1", "SGM197.2", "SGM197.3", "SGM206.1", "SGM206.2", "SGM206.3")
sintCloneDendPoints$depth = factor(sintCloneDendPoints$depth)
sintCloneDendPoints$depth = factor(sintCloneDendPoints$depth,levels(sintCloneDendPoints$depth)[c(2,1)])
sintCloneDendPoints$pop = factor(sintCloneDendPoints$pop)
sintCloneDendPoints$pop = factor(sintCloneDendPoints$pop,levels(sintCloneDendPoints$pop)[c(5, 2, 1, 4, 3)])
sintCloneDendA = ggplot() +
geom_segment(data = segment(sintCloneDData), aes(x = x, y = y, xend = xend, yend = yend2), size = 0.5) +
geom_point(data = sintCloneDendPoints, aes(x = x, y = y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
scale_fill_manual(values = gomPal, name= "Bank")+
scale_shape_manual(values = c(24, 25), name = "Depth Zone")+
geom_hline(yintercept = 0.154, color = pink, lty = 5, size = 0.75) + # creating a dashed line to indicate a clonal distance threshold
geom_text(data = subset(sintCloneDendPoints, subset = label %in% sintTechReps), aes(x = x, y = (y - .015), label = label), angle = 90) + # spacing technical replicates further from leaf
geom_text(data = subset(sintCloneDendPoints, subset = !label %in% sintTechReps), aes(x = x, y = (y - .010), label = label), angle = 90) +
labs(y = "Genetic distance (1 - IBS)") +
guides(fill = guide_legend(override.aes = list(shape = 22)))+
coord_cartesian(xlim = c(0, 220), ylim = c(0.075, 0.3), expand = FALSE) +
theme_classic()
sintCloneDend = sintCloneDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 12, color = "black", angle = 90),
axis.text.y = element_text(size = 10, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "bottom")
sintCloneDend
ggsave("../figures/extras/sintCloneDend.png", plot = sintCloneDend, height = 9, width = 30, units = "in", dpi = 300)Orbicella faveolata
ofavCloneBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "O. faveolata") %>% dplyr::filter(!Sample %in% c("OGM003", "OGM043", "OGM075", "OGM076", "OGM121"))
ofavCloneMa = as.matrix(read.table("../data/ofav/ofavClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(ofavCloneMa) = list(ofavCloneBams[,1],ofavCloneBams[,1])
ofavClonesHc = hclust(as.dist(ofavCloneMa),"ave")
ofavClonePops = ofavCloneBams$bank
ofavCloneDepth = ofavCloneBams$depthZone
ofavCloneDend = ofavCloneMa %>% as.dist() %>% hclust(.,"ave") %>% as.dendrogram()
ofavCloneDData = ofavCloneDend %>% dendro_data()
# Making the branches hang shorter so we can easily see clonal groups
ofavCloneDData$segments$yend2 = ofavCloneDData$segments$yend
for(i in 1:nrow(ofavCloneDData$segments)) {
if (ofavCloneDData$segments$yend2[i] == 0) {
ofavCloneDData$segments$yend2[i] = (ofavCloneDData$segments$y[i] - 0.01)}}
ofavCloneDendPoints = ofavCloneDData$labels
ofavCloneDendPoints$pop = ofavClonePops[order.dendrogram(ofavCloneDend)]
ofavCloneDendPoints$depth=ofavCloneDepth[order.dendrogram(ofavCloneDend)]
rownames(ofavCloneDendPoints) = ofavCloneDendPoints$label
# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(ofavCloneDData$segments)) {
if (ofavCloneDData$segments$yend[i] == 0) {
point[i] = ofavCloneDData$segments$y[i] - 0.01
} else {
point[i] = NA}}
ofavCloneDendPoints$y = point[!is.na(point)]
ofavTechReps = c("OGM024.1", "OGM024.2", "OGM024.3", "OGM066.1", "OGM066.2", "OGM066.3", "OGM108.1", "OGM108.2", "OGM108.3")
ofavCloneDendPoints$depth = factor(ofavCloneDendPoints$depth)
ofavCloneDendPoints$depth = factor(ofavCloneDendPoints$depth,levels(ofavCloneDendPoints$depth)[c(2,1)])
ofavCloneDendPoints$pop = factor(ofavCloneDendPoints$pop)
ofavCloneDendPoints$pop = factor(ofavCloneDendPoints$pop,levels(ofavCloneDendPoints$pop)[c(5, 2, 1, 4, 3)])
ofavCloneDendA = ggplot() +
geom_segment(data = segment(ofavCloneDData), aes(x = x, y = y, xend = xend, yend = yend2), size = 0.5) +
geom_point(data = ofavCloneDendPoints, aes(x = x, y = y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
scale_fill_manual(values = gomPal, name= "Bank")+
scale_shape_manual(values = c(24, 25), name = "Depth Zone")+
geom_hline(yintercept = 0.1425, color = pink, lty = 5, size = 0.75) + # creating a dashed line to indicate a clonal distance threshold
geom_text(data = subset(ofavCloneDendPoints, subset = label %in% ofavTechReps), aes(x = x, y = (y - .015), label = label), angle = 90) + # spacing technical replicates further from leaf
geom_text(data = subset(ofavCloneDendPoints, subset = !label %in% ofavTechReps), aes(x = x, y = (y - .010), label = label), angle = 90) +
labs(y = "Genetic distance (1 - IBS)") +
guides(fill = guide_legend(override.aes = list(shape = 22)))+
coord_cartesian(xlim = c(0,137), ylim = c(0.075, 0.3), expand = FALSE) +
theme_classic()
ofavCloneDend = ofavCloneDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 12, color = "black", angle = 90),
axis.text.y = element_text(size = 10, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "bottom")
ofavCloneDend
ggsave("../figures/extras/ofavCloneDend.png", plot = ofavCloneDend, height = 9, width = 30, units = "in", dpi = 300)Xestospongia muta
xestoCloneBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "X. muta") %>% dplyr::filter(!Sample %in% c("XGM013", "XGM067", "XGM129"))
xestoCloneMa = as.matrix(read.table("../data/xesto/xestoClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(xestoCloneMa) = list(xestoCloneBams[,1],xestoCloneBams[,1])
xestoClonesHc = hclust(as.dist(xestoCloneMa),"ave")
xestoClonePops = xestoCloneBams$bank
xestoCloneDepth = xestoCloneBams$depthZone
xestoCloneDend = xestoCloneMa %>% as.dist() %>% hclust(.,"ave") %>% as.dendrogram()
xestoCloneDData = xestoCloneDend %>% dendro_data()
# Making the branches hang shorter so we can easily see clonal groups
xestoCloneDData$segments$yend2 = xestoCloneDData$segments$yend
for(i in 1:nrow(xestoCloneDData$segments)) {
if (xestoCloneDData$segments$yend2[i] == 0) {
xestoCloneDData$segments$yend2[i] = (xestoCloneDData$segments$y[i] - 0.01)}}
xestoCloneDendPoints = xestoCloneDData$labels
xestoCloneDendPoints$pop = xestoClonePops[order.dendrogram(xestoCloneDend)]
xestoCloneDendPoints$depth=xestoCloneDepth[order.dendrogram(xestoCloneDend)]
rownames(xestoCloneDendPoints) = xestoCloneDendPoints$label
# Making points at the leaves to place symbols for populations
point = as.vector(NA)
for(i in 1:nrow(xestoCloneDData$segments)) {
if (xestoCloneDData$segments$yend[i] == 0) {
point[i] = xestoCloneDData$segments$y[i] - 0.01
} else {
point[i] = NA}}
xestoCloneDendPoints$y = point[!is.na(point)]
xestoTechReps = c("XGM034.1", "XGM034.2", "XGM034.3", "XGM072.1", "XGM072.2", "XGM072.3", "XGM158.1", "XGM158.2", "XGM158.3", "XGM179.1", "XGM179.2", "XGM179.3", "XGM193.1", "XGM193.2", "XGM193.3", "XGM203.1", "XGM203.2", "XGM203.3")
xestoCloneDendPoints$depth = factor(xestoCloneDendPoints$depth)
xestoCloneDendPoints$depth = factor(xestoCloneDendPoints$depth,levels(xestoCloneDendPoints$depth)[c(2,1)])
xestoCloneDendPoints$pop = factor(xestoCloneDendPoints$pop)
xestoCloneDendPoints$pop = factor(xestoCloneDendPoints$pop,levels(xestoCloneDendPoints$pop)[c(5, 2, 1, 4, 3)])
xestoCloneDendA = ggplot() +
geom_segment(data = segment(xestoCloneDData), aes(x = x, y = y, xend = xend, yend = yend2), size = 0.5) +
geom_point(data = xestoCloneDendPoints, aes(x = x, y = y, fill = pop, shape = depth), size = 4, stroke = 0.25) +
scale_fill_manual(values = gomPal, name= "Bank")+
scale_shape_manual(values = c(24, 25), name = "Depth Zone")+
geom_hline(yintercept = 0.174, color = pink, lty = 5, linewidth = 0.75) + # creating a dashed line to indicate a clonal distance threshold
geom_text(data = subset(xestoCloneDendPoints, subset = label %in% xestoTechReps), aes(x = x, y = (y - .015), label = label), angle = 90) + # spacing technical replicates further from leaf
geom_text(data = subset(xestoCloneDendPoints, subset = !label %in% xestoTechReps), aes(x = x, y = (y - .010), label = label), angle = 90) +
labs(y = "Genetic distance (1 - IBS)") +
guides(fill = guide_legend(override.aes = list(shape = 22)))+
coord_cartesian(xlim = c(0, 217), ylim = c(0.075, 0.3), expand = FALSE) +
theme_classic()
xestoCloneDend = xestoCloneDendA + theme(
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(size = 12, color = "black", angle = 90),
axis.text.y = element_text(size = 10, color = "black"),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
panel.grid = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
plot.background = element_blank(),
legend.key = element_blank(),
legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.position = "bottom")
# xestoCloneDend
ggsave("../figures/extras/xestoCloneDend.png", plot = xestoCloneDend, height = 8, width = 35, units = "in", dpi = 300)Montastraea cavernosa
Dendrogram and structure
mcavBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "M. cavernosa") %>% dplyr::filter(!Sample %in% c("MGM072.1", "MGM072.3", "MGM013.2", "MGM013.3", "MGM008.1", "MGM008.2", "MGM010", "MGM024.1", "MGM004", "MGM002.2", "MGM002.3", "MGM038.2", "MGM038.3"))# list of bams files and their populations (tech reps removed)
mcavBams$Sample <- gsub("\\.[1-3]*$", "", mcavBams$Sample)
mcavBams = mcavBams %>% dplyr::mutate(sample = Sample) %>% dplyr::select(-Sample) %>% dplyr::relocate(sample)
mcavBams$bank = factor(mcavBams$bank)
mcavBams$bank = factor(mcavBams$bank, levels = levels(mcavBams$bank)[c(5, 2, 1, 3, 4)])
mcavBams$depthZone = factor(mcavBams$depthZone)
mcavBams$depthZone = factor(mcavBams$depthZone, levels = levels(mcavBams$depthZone)[c(2, 1)])
mcavMa = as.matrix(read.table("../data/mcav/mcavNoClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(mcavMa) = list(mcavBams[,3],mcavBams[,3])
## admixture K = 2
#-------------------------------------
mcavK2ad = read.table("../data/mcav/admix/mcavK2.output") %>% dplyr::select(V6, V7)
mcavK2ad %>% summarise(sum(V6),sum(V7))## sum(V6) sum(V7)
## 1 119.4213 49.5787mcavAdmixK2 = mcavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(mcavK2ad) %>%
dplyr::rename("Mc1" = "V7", "Mc2" = "V6") %>% dplyr::select(order(colnames(.))) %>% dplyr::relocate(sample)
mcavAdmixK2_melt = melt(mcavAdmixK2, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 3
#-------------------------------------
mcavK3ad = read.table("../data/mcav/admix/mcavK3.output") %>% dplyr::select(V6, V7, V8)
mcavK3ad %>% summarise(sum(V6),sum(V7), sum(V8))## sum(V6) sum(V7) sum(V8)
## 1 65.4637 43.7152 59.8217mcavAdmixK3 = mcavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(mcavK3ad) %>%
dplyr::rename("Mc1" = "V7", "Mc2" = "V6",, "Mc3" = "V8") %>%
dplyr::select(order(colnames(.)))
mcavAdmixK3_melt = melt(mcavAdmixK3, id = c("sample", "bank", "depth", "depthm"))
## ngsAdmix K = 4
#-------------------------------------
mcavK4ad = read.table("../data/mcav/admix/mcavK4.output") %>% dplyr::select(V6, V7, V8, V9)
mcavK4ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9))## sum(V6) sum(V7) sum(V8) sum(V9)
## 1 43.3693 39.3483 41.005 45.2775mcavAdmixK4 = mcavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(mcavK4ad) %>%
dplyr::rename("Mc1" = "V7", "Mc2" = "V6", "Mc3" = "V8", "Mc4" = "V9") %>%
dplyr::select(order(colnames(.)))
mcavAdmixK4_melt = melt(mcavAdmixK4, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 5
#-------------------------------------
mcavK5ad = read.table("../data/mcav/admix/mcavK5.output") %>% dplyr::select(V6, V7, V8, V9, V10)
mcavK5ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10)
## 1 39.392 33.285 36.313 35.6101 24.3984mcavAdmixK5 = mcavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(mcavK5ad) %>%
dplyr::rename("Mc1" = "V7", "Mc2" = "V6", "Mc3" = "V8", "Mc4" = "V9", "Mc5" = "V10") %>% dplyr::select(order(colnames(.)))
mcavAdmixK5_melt = melt(mcavAdmixK5, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 6
#-------------------------------------
mcavK6ad = read.table("../data/mcav/admix/mcavK6.output") %>% dplyr::select(V6, V7, V8, V9, V10,V11)
mcavK6ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10), sum(V11))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10) sum(V11)
## 1 32.9952 29.8928 30.4508 29.7122 20.8536 25.0957mcavAdmixK6 = mcavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(mcavK6ad) %>%
dplyr::rename("Mc1" = "V7", "Mc2" = "V6", "Mc3" = "V8", "Mc4" = "V9", "Mc5" = "V10", "Mc6" = "V11") %>% dplyr::select(order(colnames(.)))
mcavAdmixK6_melt = melt(mcavAdmixK6, id = c("sample", "bank", "depth", "depthm"))
## construct figure
#-------------------------------------
# ggtree(tr) + geom_nodelab(aes(label = node), hjust = -0.5)
{
mcavTr = hclust(dist(mcavMa),"ave")
mcavP1 = ggtree(mcavTr, layout = "rectangular", size = 0.35)
mcavP2 = facet_plot(mcavP1, panel = "location", data=mcavAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = bank), color = "grey25", size = 0.1) +
scale_fill_manual("Site", values = gomPal, guide = guide_legend(order = 1)) +
new_scale_fill()
mcavP3 = facet_plot(mcavP2, panel = "depth zone", data=mcavAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depth), color = "grey25", size = 0.1) +
scale_fill_manual("Depth zone", values = c("#A7E1BCFF", "#414083FF"), guide = guide_legend(order = 2)) +
new_scale_fill()
mcavP4 = facet_plot(mcavP3, panel = "depth", data=mcavAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depthm), color = "grey25", size = 0.1) +
scale_fill_viridis_c("Depth (m)", option = "mako", trans = "reverse", limits = c(60, 10)) +
new_scale_fill()
mcavP5 = facet_plot(mcavP4, panel="K = 2", data=mcavAdmixK2_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
scale_fill_manual("Lineage",values = mcavKColPal[1:6])
mcavP6 = facet_plot( mcavP5, panel="K = 3", data=mcavAdmixK3_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
mcavP7 = facet_plot(mcavP6, panel="K = 4", data=mcavAdmixK4_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
mcavP8 = facet_plot(mcavP7, panel="K = 5", data=mcavAdmixK5_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
mcavP9 = facet_plot(mcavP8, panel="K = 6", data=mcavAdmixK6_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
theme_tree(strip.text = element_blank(),
panel.spacing = unit(.1, "line"))
}
mcavImg = image_read("../figures/icons/mcav.png") %>% image_ggplot()
mcavStructure = facet_widths(mcavP9, widths = c(0.5, 0.025, 0.025, 0.025, 0.2, 0.1, 0.1, 0.1, 0.1)) #+ inset_element(mcavImg, 0.0, 0.9, 0.1, 1.0, align_to = "full", ignore_tag = TRUE)
mcavStructure
Population structure
mcavData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "M. cavernosa") %>% dplyr::filter(!Sample %in% c("MGM072.1", "MGM072.3", "MGM013.2", "MGM013.3", "MGM008.1", "MGM008.2", "MGM010", "MGM024.1", "MGM004", "MGM002.2", "MGM002.3", "MGM038.2", "MGM038.3"))
mcavData$depthZone = factor(mcavData$depthZone)
mcavData$depthZone = factor(mcavData$depthZone, levels(mcavData$depthZone)[c(2,1)])
mcavData$bank = factor(mcavData$bank)
mcavData$bank = factor(mcavData$bank,levels(mcavData$bank)[c(5, 2, 1, 4, 3)])
mcavPcadmix = read.table("../data/mcav/admix/mcavK2.output") %>%dplyr::select(V6, V7)
mcavPcadmix %>% summarise(sum(V6),sum(V7)) ## sum(V6) sum(V7)
## 1 119.4213 49.5787mcavPcadmix = mcavData %>% cbind(mcavPcadmix) %>% rename("cluster1" = "V7", "cluster2" = "V6") %>%dplyr::select(order(colnames(.)))
mcavPcadmixClust = mcavPcadmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, 0))))
mcavPcadmix = mcavPcadmix %>% mutate(mcavPcadmixClust)
mcavPcadmix$cluster = as.factor(mcavPcadmix$cluster)
levels(mcavPcadmix$cluster) = c("Mc_1", "Mc_2", "Admixed")
mcavSiteLineages = mcavPcadmix %>% dplyr::select(depthZone, cluster) %>%
group_by(depthZone) %>% count(cluster) %>% mutate(Freq = n/sum(n)) %>% apply(2, function(x) gsub("\\s+", "", x)) %>% as.data.frame()PCA
mcavCov = read.table("../data/mcav/mcavPcangsd.cov") %>% as.matrix()
mcavPcAdmix = read.table("../data/mcav/admix/mcavK2.output") %>%dplyr::select(V6, V7)
mcavPcAdmix %>% summarise(sum(V6),sum(V7)) ## sum(V6) sum(V7)
## 1 119.4213 49.5787mcavPcAdmix = mcavPcAdmix %>% rename("cluster1" = "V7", "cluster2" = "V6") %>% dplyr::select(order(colnames(.)))
mcavPcaEig = eigen(mcavCov)
mcavPcaVar = mcavPcaEig$values/sum(mcavPcaEig$values)*100
head(mcavPcaVar)## [1] 2.9308183 0.8264292 0.8100717 0.8070106 0.7996972 0.7900787mcavPcangsd = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "M. cavernosa") %>% dplyr::filter(!Sample %in% c("MGM072.1", "MGM072.3", "MGM013.2", "MGM013.3", "MGM008.1", "MGM008.2", "MGM010", "MGM024.1", "MGM004", "MGM002.2", "MGM002.3", "MGM038.2", "MGM038.3")) %>% dplyr::select("sample" = Sample, "bank" = bank, "depth" = depthZone, "depthm" = depthM)
mcavPcangsd$sitedepth = as.factor(paste(mcavPcangsd$bank, mcavPcangsd$depth, sep = " "))
mcavPcangsd$sitedepth = factor(mcavPcangsd$sitedepth, levels(mcavPcangsd$sitedepth)[c(7, 6, 3, 2, 1, 4, 5)])
mcavPcangsd$bank = factor(mcavPcangsd$bank)
mcavPcangsd$bank = factor(mcavPcangsd$bank, levels(mcavPcangsd$bank)[c(5, 2, 1, 3, 4)])
mcavPcangsd$depth = factor(mcavPcangsd$depth)
mcavPcangsd$depth = factor(mcavPcangsd$depth, levels(mcavPcangsd$depth)[c(2, 1)])
mcavPcangsd$PC1 = mcavPcaEig$vectors[,1]
mcavPcangsd$PC2 = mcavPcaEig$vectors[,2]
mcavPcangsd$PC3 = mcavPcaEig$vectors[,3]
mcavPcangsd$PC4 = mcavPcaEig$vectors[,4]
mcavPcangsdClust = mcavPcAdmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, 0))))
# pcangsdClust$clusterX = as.vector(d_clust$classification)
mcavPcangsd = mcavPcangsd %>% mutate(mcavPcangsdClust)
mcavPcangsd$cluster = as.factor(mcavPcangsd$cluster)
levels(mcavPcangsd$cluster) = c("Mc_Deep", "Mc_1", "Admixed")
mcavBamsClusters = mcavPcangsd %>% dplyr::select(sample, cluster) %>% dplyr::arrange(sample)
mcavBamsSamples = read.delim("../data/mcav/bamsNoClones", header = FALSE)
mcavBamsClusters$sample = mcavBamsSamples$V1
# bamsClusters = as.data.frame(bamsClusters)
write.table(x = mcavBamsClusters, file = "../data/mcav/mcavBamsClusters", sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
mcavPcangsd = merge(mcavPcangsd, aggregate(cbind(mean.x = PC1, mean.y = PC2)~sitedepth, mcavPcangsd, mean), by="sitedepth")mcavPcangsd %>% dplyr::group_by(depth,cluster) %>% dplyr::summarize(n = n())## `summarise()` has grouped output by 'depth'. You can override using the `.groups`
## argument.## # A tibble: 5 × 3
## # Groups: depth [2]
## depth cluster n
## <fct> <fct> <int>
## 1 Shallow Mc_1 45
## 2 Shallow Admixed 8
## 3 Mesophotic Mc_Deep 14
## 4 Mesophotic Mc_1 57
## 5 Mesophotic Admixed 45Plot PCA
mcavPcaPlot12SA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.25) +
geom_vline(xintercept = 0, color = "gray90", size = 0.25) +
geom_point(data = mcavPcangsd, aes(x = PC1, y = PC2, fill = bank, shape = depth, color = bank), stroke = 0, size = 3, alpha = 0.5, show.legend = FALSE) +
geom_point(data = mcavPcangsd, aes(x = mean.x, y = mean.y, fill = bank, shape = depth), color = "black", size = 3.75, alpha = 1, stroke = 0.25) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = gomPal, name = "Site") +
scale_color_manual(values = gomPal, name = "Site") +
labs(x = paste0("PC 1 (", format(round(mcavPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(mcavPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = guide_legend(override.aes = list(size = 2, stroke = 0.25, alpha = NA), order = 2, ncol = 1), fill = guide_legend(override.aes = list(shape = 22, size = 2, fill = gomPal[c(1,2,3,4,5)], alpha = NA), order = 1, ncol = 1)) +
theme_bw()
mcavPcaPlot12S = mcavPcaPlot12SA +
pcaTheme +
theme(legend.position = c(0.87, 0.76))
mcavPcaPlot12S
mcavPcaPlot12LA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
geom_point(data = mcavPcangsd, aes(x = PC1, y = PC2, fill = cluster, shape = depth), color = "black", size = 3, alpha = 1, show.legend = TRUE) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = mcavKColPal[c(1, 2, 7)], name = "Lineage") +
labs(x = paste0("PC 1 (", format(round(mcavPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(mcavPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = "none", fill = guide_legend(override.aes = list(shape = 22, size = 2, fill = mcavKColPal[c(1, 2, 7)], alpha = NA), order = 1, ncol = 1))+
theme_bw()
mcavPcaPlot12L = mcavPcaPlot12LA +
pcaTheme +
theme(legend.position = c(0.9,0.86))
mcavPcaPlot12L
Put all plots together
mcavPcaPlots = ((mcavPcaPlot12S + theme(axis.title.y = element_text(margin = ggplot2::margin(r = -20, unit = "pt")))) | mcavPcaPlot12L) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
panel.background = element_rect(fill = "white"),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
mcavPcaPlots
Admixture
Prepare admixture outputs for plotting
mcavAdmix = mcavPcangsd %>%dplyr::select(-PC1, -PC2, -PC3, -PC4, -cluster, -depthm, -mean.x, -mean.y)
mcavAdmix$bank = factor(mcavAdmix$bank)
mcavAdmix = arrange(mcavAdmix, bank, depth, -cluster1, -cluster2)
mcavPopCounts = mcavAdmix %>% group_by(bank, depth) %>% summarise(n = n())## `summarise()` has grouped output by 'bank'. You can override using the `.groups`
## argument.mcavPopCounts## # A tibble: 7 × 3
## # Groups: bank [5]
## bank depth n
## <fct> <fct> <int>
## 1 WFGB Shallow 27
## 2 WFGB Mesophotic 27
## 3 EFGB Shallow 26
## 4 EFGB Mesophotic 21
## 5 Bright Mesophotic 28
## 6 Geyer Mesophotic 12
## 7 McGrail Mesophotic 28mcavPopCountList = reshape2::melt(lapply(mcavPopCounts$n,function(x){c(1:x)}))
mcavAdmix$ord = mcavPopCountList$value
mcavAdmixMelt = melt(mcavAdmix, id.vars=c("sample", "bank", "depth", "sitedepth", "ord"), variable.name="Ancestry", value.name="Fraction")
mcavAdmixMelt$Ancestry = factor(mcavAdmixMelt$Ancestry)
mcavAdmixMelt$Ancestry = factor(mcavAdmixMelt$Ancestry, levels = rev(levels(mcavAdmixMelt$Ancestry)))
mcavPopAnno = data.frame(x1 = c(0.5, 0.5, 0.5, 0.5, 0.5), x2 = c(28.5, 28.5, 28.5, 28.5, 28.5),
y1 = -0.1, y2 = -0.1, sample = NA, Ancestry = NA, depth = "Mesophotic",
ord = NA, Fraction = NA,
bank = c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"))
mcavPopAnno$bank = factor(mcavPopAnno$bank)
mcavPopAnno$bank = factor(mcavPopAnno$bank, levels = levels(mcavPopAnno$bank)[c(5, 2, 1, 3, 4)])Make admixture plots
mcavAdmixPlotA = ggplot(data = mcavAdmixMelt, aes(x = ord, y = Fraction, fill = Ancestry, order = sample)) +
geom_segment(data = mcavPopAnno, aes(x = x1, xend = x2, y = -.1, yend = -.1, color = bank), size = 7) +
geom_bar(stat = "identity", position = "fill", width = 1, colour = "grey25", size = 0.2) +
facet_grid(factor(depth) ~ bank, switch = "both", scales = "free_x") +
geom_text(data = (mcavAdmixMelt %>% filter(depth == "Mesophotic", bank %in% c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"), sample %in% c(
"MGM116", "MGM015", "MGM001", "MGM144", "MGM062"), Ancestry == "cluster1")), x = 15.5, y = -.1, aes(label = bank), size = 4, color = "black") +
scale_fill_manual(values = mcavKColPal[1:2]) +
scale_color_manual(values = gomPal) +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
coord_cartesian(ylim = c(0.0, 1.0), clip = "off") +
theme_bw()
mcavAdmixPlot = mcavAdmixPlotA +
theme_bw()+
admixTheme
mcavAdmixPlot
Lineage demographics
Lineage depth distribuiton
leveneTest(lm(depthm ~ cluster, data = mcavPcangsd))## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 20.177 0.00000001432 ***
## 166
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mcavDepthAov = welch_anova_test(depthm ~ cluster, data = subset(mcavPcangsd, subset = mcavPcangsd$cluster!="Admixed"))
mcavDF = mcavDepthAov$statistic[[1]]
mcavDepthPH = games_howell_test(depthm ~ cluster, data = mcavPcangsd, conf.level = 0.95) %>% as.data.frame()
mcavLineageViolinA = ggplot(data = mcavPcangsd, aes(x = cluster, y = depthm)) +
annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 30, ymax = Inf, fill = "black", alpha = 0.10, color = NA) +
geom_violin(aes(fill = cluster),adjust = 1, linewidth = 0, color = "black", alpha = 0.75, width = 0.9, trim = F, scale = "area") +
geom_beeswarm(aes(fill = cluster), shape = 21, size = 2, cex = 1.5, alpha = 0.5) +
geom_violin(adjust = 1, linewidth = 0.4, color = "black", alpha = 1, width = 0.9, trim = F, fill = NA, scale = "area") +
geom_boxplot(width = 0.1, color = "black", fill = "white", outlier.colour = NA, linewidth = 0.6, alpha = 0.5) +
annotate(geom = "text", x = 2.75, y =65, label = bquote(italic("F")~" = "~.(mcavDF)*"; "~italic("p")~" < 0.001"), size = 3) +
scale_fill_discrete(type = mcavKColPal[c(1, 2, 7)], name = "Lineage") +
scale_color_discrete(type = mcavKColPal[c(1, 2, 7)], name = "Lineage") +
xlab("Lineage") +
ylab("Depth (m)") +
scale_y_reverse(breaks = seq(5, 90, 5)) +
theme_bw()
mcavLineageViolin = mcavLineageViolinA + violinTheme
mcavLineageViolin
Symbiodiniaceae
mcavData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "M. cavernosa") %>% dplyr::filter(!Sample %in% c("MGM072.1", "MGM072.3", "MGM013.2", "MGM013.3", "MGM008.1", "MGM008.2", "MGM010", "MGM024.1", "MGM004", "MGM002.2", "MGM002.3", "MGM038.2", "MGM038.3")) %>% dplyr::select("sample" = Sample, "site" = bank, "depth" = depthZone)
mcavZoox = read.delim("../data/mcav/mcavZooxReads", header = FALSE, check.names = FALSE)
head(mcavZoox)## V1 V2
## 1 MGM001.trim.zoox.bt2.bam NA
## 2 chr1 15
## 3 chr2 32
## 4 chr3 12
## 5 chr4 29
## 6 chr5 113# Reconstruct read mapping output into dataframe usable for analysis
mcavZoox$V2[is.na(mcavZoox$V2)] <- as.character(mcavZoox$V1[is.na(mcavZoox$V2)])
mcavZoox$V1 = gsub(pattern = "MGM*.", "chr", mcavZoox$V1)
mcavZoox$V2 = gsub(".trim.*", "", mcavZoox$V2)
mcavZoox = mcavZoox %>% filter(mcavZoox$V1 != "*")
mcavZooxLst = split(mcavZoox$V2, as.integer(gl(length(mcavZoox$V2), 20, length(mcavZoox$V2))))
mcavZooxMaps = NULL
for(i in mcavZooxLst){
mcavZooxMaps = rbind(mcavZooxMaps, data.frame(t(i)))
}
# remove tech reps
mcavZooxMaps = mcavZooxMaps %>% dplyr::filter(!X1 %in% c("MGM072.1", "MGM072.3", "MGM013.2", "MGM013.3", "MGM008.1", "MGM008.2", "MGM010", "MGM024.1", "MGM004", "MGM002.2", "MGM002.3", "MGM038.2", "MGM038.3"))
colnames(mcavZooxMaps) = c("sample",mcavZoox$V1[c(2:20)])
# convert characters to numeric
str(mcavZooxMaps)## 'data.frame': 169 obs. of 20 variables:
## $ sample: chr "MGM001" "MGM002.1" "MGM003" "MGM005" ...
## $ chr1 : chr "15" "88" "38" "108" ...
## $ chr2 : chr "32" "98" "69" "125" ...
## $ chr3 : chr "12" "83" "46" "79" ...
## $ chr4 : chr "29" "94" "36" "101" ...
## $ chr5 : chr "113" "133" "44" "148" ...
## $ chr6 : chr "51" "108" "44" "99" ...
## $ chr7 : chr "24" "31" "41" "72" ...
## $ chr8 : chr "77" "87" "87" "141" ...
## $ chr9 : chr "55" "74" "66" "141" ...
## $ chr10 : chr "718" "2764" "1330" "6699" ...
## $ chr11 : chr "881" "3522" "1501" "8839" ...
## $ chr12 : chr "1065" "4409" "1835" "9631" ...
## $ chr13 : chr "951" "3580" "1579" "8742" ...
## $ chr14 : chr "1038" "3669" "1738" "8645" ...
## $ chr15 : chr "238" "642" "382" "1639" ...
## $ chr16 : chr "26" "60" "96" "238" ...
## $ chr17 : chr "29" "69" "40" "81" ...
## $ chr18 : chr "13" "28" "27" "90" ...
## $ chr19 : chr "5" "4" "6" "21" ...for(i in c(2:20)){
mcavZooxMaps[,i] = as.numeric(mcavZooxMaps[,i])
}
str(mcavZooxMaps)## 'data.frame': 169 obs. of 20 variables:
## $ sample: chr "MGM001" "MGM002.1" "MGM003" "MGM005" ...
## $ chr1 : num 15 88 38 108 56 88 34 22 56 11 ...
## $ chr2 : num 32 98 69 125 69 41 36 34 65 12 ...
## $ chr3 : num 12 83 46 79 44 33 28 21 72 15 ...
## $ chr4 : num 29 94 36 101 66 45 30 26 65 13 ...
## $ chr5 : num 113 133 44 148 122 116 128 3 9 3 ...
## $ chr6 : num 51 108 44 99 72 72 25 9 8 5 ...
## $ chr7 : num 24 31 41 72 36 41 29 5 1 3 ...
## $ chr8 : num 77 87 87 141 82 58 52 9 13 6 ...
## $ chr9 : num 55 74 66 141 100 67 48 15 14 17 ...
## $ chr10 : num 718 2764 1330 6699 4004 ...
## $ chr11 : num 881 3522 1501 8839 3526 ...
## $ chr12 : num 1065 4409 1835 9631 4135 ...
## $ chr13 : num 951 3580 1579 8742 3537 ...
## $ chr14 : num 1038 3669 1738 8645 3452 ...
## $ chr15 : num 238 642 382 1639 726 ...
## $ chr16 : num 26 60 96 238 331 68 10 6 13 15 ...
## $ chr17 : num 29 69 40 81 144 64 9 2 9 8 ...
## $ chr18 : num 13 28 27 90 231 31 5 7 11 12 ...
## $ chr19 : num 5 4 6 21 29 2 1 2 1 2 ...# collapse fake chromosomes into representative genera
mcavZooxMaps$Symbiodinium = rowSums(mcavZooxMaps[2:6])
mcavZooxMaps$Breviolum = rowSums(mcavZooxMaps[7:10])
mcavZooxMaps$Cladocopium = rowSums(mcavZooxMaps[11:16])
mcavZooxMaps$Durusdinium = rowSums(mcavZooxMaps[17:20])
# keep genera totals and turn into proportions for barplot
mcavZooxMaps = mcavZooxMaps[,c(1, 21:24)]
mcavZooxProp = mcavZooxMaps
mcavZooxProp$sum = apply(mcavZooxProp[, c(2:length(mcavZooxProp[1,]))], 1, function(x) {
sum(x, na.rm = T)
})
mcavZooxProp = cbind(mcavZooxProp$sample, (mcavZooxProp[, c(2:(ncol(mcavZooxProp)-1))]
/ mcavZooxProp$sum))
colnames(mcavZooxProp)[1] = "sample"
head(mcavZooxProp)## sample Symbiodinium Breviolum Cladocopium Durusdinium
## 1 MGM001 0.037416232 0.038533135 0.9104617 0.013588980
## 2 MGM002.1 0.025379931 0.015350765 0.9510311 0.008238244
## 3 MGM003 0.025874514 0.026429761 0.9289284 0.018767351
## 4 MGM005 0.012292119 0.009925721 0.9683604 0.009421766
## 5 MGM006 0.017194875 0.013967826 0.9334361 0.035401214
## 6 MGM007 0.006753366 0.004976165 0.9848206 0.003449862# Check that all samples total to 1
apply(mcavZooxProp[, c(2:(ncol(mcavZooxProp)))], 1, function(x) {
sum(x, na.rm = T)
})## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [42] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [124] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [165] 1 1 1 1 1# add sample metadata to proportions
mcavSym = mcavData %>% left_join(mcavZooxProp)## Joining with `by = join_by(sample)`mcavAdmixOrd = mcavAdmix %>% dplyr::select(sample, ord)
mcavPcangsdITS =mcavPcangsd %>% dplyr::select(sample, cluster)
mcavSym = mcavSym %>% left_join(mcavAdmixOrd) %>% relocate(ord,.after = depth) %>% left_join(mcavPcangsdITS) %>% relocate(cluster, .after = depth)## Joining with `by = join_by(sample)`
## Joining with `by = join_by(sample)`set.seed(694)
mcavSymPerm = adonis2(decostand(mcavSym[, c(6:((ncol(mcavSym))))], "normalize") ~ site * depth * cluster, data = mcavSym, permutations = 9999, method = "bray")
mcavSymPerm## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = decostand(mcavSym[, c(6:((ncol(mcavSym))))], "normalize") ~ site * depth * cluster, data = mcavSym, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Model 17 0.049142 0.20293 2.2614 0.0183 *
## Residual 151 0.193022 0.79707
## Total 168 0.242164 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mcavSym$depth = factor(mcavSym$depth)
mcavSym$depth = factor(mcavSym$depth, levels = levels(mcavSym$depth)[c(2, 1)])
mcavSym$site = factor(mcavSym$site)
mcavSym$site = factor(mcavSym$site, levels = levels(mcavSym$site)[c(5, 2, 1, 3, 4)])
#turn into melted dataframe with otustack() and remove "summ" rows
mcavGssSym = otuStack(mcavSym, count.columns = c(6:length(mcavSym[1, ])),
condition.columns = c(1:5)) %>% filter(otu != "summ") %>% droplevels()
#check that levels are correct/ordered
levels(mcavGssSym$otu)## [1] "Symbiodinium" "Breviolum" "Cladocopium" "Durusdinium"levels(mcavGssSym$depth)## [1] "Shallow" "Mesophotic"levels(mcavGssSym$site)## [1] "WFGB" "EFGB" "Bright" "Geyer" "McGrail"Symbiodiniaceae genera barplot
mcavZooxAnno = data.frame(x1 = c(0.5, 0.5, 0.5, 0.5, 0.5), x2 = c(28.5, 28.5, 28.5, 28.5, 28.5),
y1 = -0.22, y2 = -0.22, site = c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"))
mcavZooxAnno$site = factor(mcavZooxAnno$site)
mcavZooxAnno$site = factor(mcavZooxAnno$site, levels = levels(mcavZooxAnno$site)[c(5, 2, 1, 3, 4)])
mcavGssSymPlot = mcavGssSym %>% left_join(mcavZooxAnno, by = "site") %>% left_join(mcavPcangsdITS)## Joining with `by = join_by(sample, cluster)`mcavGssSymPlot$ord = as.numeric(mcavGssSymPlot$ord)
mcavZooxPlotA = ggplot(data = mcavGssSymPlot, aes(x = ord, y = count, fill = otu, order = ord)) +
geom_point(aes(x=1, y=0.5, fill = otu), shape = 22, size = 0) +
geom_bar(stat = "identity", position = "stack", colour = "grey25", width = 1, size = 0.2, show.legend = FALSE) +
xlab("Population") +
scale_fill_manual(values = colPalZoox, name = "Symbiodiniaceae \ngenus") +
geom_segment(data = (mcavGssSymPlot %>% filter(depth == "Mesophotic")), aes(x = x1, xend = x2, y = y1, yend = y2, color = site), size = 7) +
scale_color_manual(values = gomPal, guide = "none") +
ggnewscale::new_scale_color() +
geom_segment(data = mcavGssSymPlot, aes(x = ord-0.5, xend = ord+0.5, color = cluster), y = -.07, yend = -.07, linewidth = 4) +
scale_color_manual(values = mcavKColPal[c(1, 2, 7)], name = "Lineage") +
coord_cartesian(ylim = c(0, 1), xlim = c(0.5, 28.5), clip = "off") +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
facet_grid(factor(depth) ~ site, drop = TRUE, scales = "free", switch = "both", space = "free") +
geom_text(data = subset(mcavGssSymPlot, subset = otu == "Symbiodinium") %>% filter(sample %in% c("MGM116", "MGM015", "MGM001", "MGM144", "MGM062")), x = 15.5, y = -.205, aes(label = site), size = 4.4, color = "#000000") +
guides(fill = guide_legend(override.aes = list(size = 4), ncol = 1, order = 1), color = guide_legend(override.aes = list(size = .1), ncol = 1)) +
theme_bw()
mcavZooxPlot = mcavZooxPlotA + zooxTheme
mcavZooxPlot
All M. cavernosa plots
mcavPlot1 = mcavStructure + inset_element(mcavImg, -.1,.7,.3,1, align_to = "full", ignore_tag = TRUE)
mcavPlot2 = (mcavPcaPlots + mcavLineageViolin) + plot_layout(guides = "collect")
mcavPlot3 = mcavAdmixPlot + mcavZooxPlot
mcavPlots = mcavPlot1 + mcavPlot2 + mcavPlot3 +
plot_layout(heights = c(2,1,1)) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
ggsave("../figures/figure2.png", plot = mcavPlots, height = 11, width = 12, units = "in", dpi = 300)Stephanocoenia intersepta
Dendrogram and structure
sintBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "S. intersepta") %>% dplyr::filter(!Sample %in% c("SGM007", "SGM010", "SGM013", "SGM027.2", "SGM027.3", "SGM046.1", "SGM046.3", "SGM063", "SGM152.2", "SGM152.3", "SGM169", "SGM177", "SGM179", "SGM186.1", "SGM186.3", "SGM197.2", "SGM197.3", "SGM206.2", "SGM206.3"))
sintBams$Sample <- gsub("\\.[1-3]*$", "", sintBams$Sample)
sintBams = sintBams %>% dplyr::mutate(sample = Sample) %>% dplyr::select(-Sample) %>% dplyr::relocate(sample)
sintBams$bank = factor(sintBams$bank)
sintBams$bank = factor(sintBams$bank, levels = levels(sintBams$bank)[c(5, 2, 1, 3, 4)])
sintBams$depthZone = factor(sintBams$depthZone)
sintBams$depthZone = factor(sintBams$depthZone, levels = levels(sintBams$depthZone)[c(2, 1)])
sintMa = as.matrix(read.table("../data/sint/sintNoClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(sintMa) = list(sintBams[,3],sintBams[,3])
## admixture K = 2
#-------------------------------------
sintK2ad = read.table("../data/sint/admix/sintK2.output") %>% dplyr::select(V6, V7)
sintK2ad %>% summarise(sum(V6),sum(V7))## sum(V6) sum(V7)
## 1 103.7485 99.2515sintAdmixK2 = sintBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(sintK2ad) %>%
dplyr::rename("Si1" = "V6", "Si2" = "V7")
sintAdmixK2_melt = melt(sintAdmixK2, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 3
#-------------------------------------
sintK3ad = read.table("../data/sint/admix/sintK3.output") %>% dplyr::select(V6, V7, V8)
sintK3ad %>% summarise(sum(V6),sum(V7), sum(V8))## sum(V6) sum(V7) sum(V8)
## 1 73.0318 70.4823 59.4852sintAdmixK3 = sintBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(sintK3ad) %>%
dplyr::rename("Si1" = "V6", "Si2" = "V7", "Si3" = "V8") %>%
dplyr::select(order(colnames(.)))
sintAdmixK3_melt = melt(sintAdmixK3, id = c("sample", "bank", "depth", "depthm"))
## ngsAdmix K = 4
#-------------------------------------
sintK4ad = read.table("../data/sint/admix/sintK4.output") %>% dplyr::select(V6, V7, V8, V9)
sintK4ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9))## sum(V6) sum(V7) sum(V8) sum(V9)
## 1 50.0275 48.5974 48.3334 56.0415sintAdmixK4 = sintBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(sintK4ad) %>%
dplyr::rename("Si1" = "V6", "Si2" = "V7", "Si3" = "V8", "Si4" = "V9") %>%
dplyr::select(order(colnames(.)))
sintAdmixK4_melt = melt(sintAdmixK4, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 5
#-------------------------------------
sintK5ad = read.table("../data/sint/admix/sintK5.output") %>% dplyr::select(V6, V7, V8, V9, V10)
sintK5ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10)
## 1 43.6614 41.4828 40.1593 44.0573 33.6384sintAdmixK5 = sintBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(sintK5ad) %>%
dplyr::rename("Si1" = "V6", "Si2" = "V7", "Si3" = "V8", "Si4" = "V9", "Si5" = "V10") #%>%
sintAdmixK5_melt = melt(sintAdmixK5, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 6
#-------------------------------------
sintK6ad = read.table("../data/sint/admix/sintK6.output") %>% dplyr::select(V6, V7, V8, V9, V10,V11)
sintK6ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10), sum(V11))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10) sum(V11)
## 1 30.0317 36.7717 36.8505 36.0188 32.8497 30.4772sintAdmixK6 = sintBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(sintK6ad) %>%
dplyr::rename("Si1" = "V6", "Si2" = "V7", "Si3" = "V8", "Si4" = "V9", "Si5" = "V10", "Si6" = "V11") #%>%
sintAdmixK6_melt = melt(sintAdmixK6, id = c("sample", "bank", "depth", "depthm"))
## construct figure
#-------------------------------------
# ggtree(tr) + geom_nodelab(aes(label = node), hjust = -0.5)
{
sintTr = hclust(dist(sintMa),"ave")
sintP1 = ggtree(sintTr, layout = "rectangular", size = 0.35)
sintP2 = facet_plot(sintP1, panel = "location", data=sintAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = bank), color = "grey25", size = 0.1) +
scale_fill_manual("Site", values = gomPal, guide = guide_legend(order = 1)) +
new_scale_fill()
sintP3 = facet_plot(sintP2, panel = "depth zone", data=sintAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depth), color = "grey25", size = 0.1) +
scale_fill_manual("Depth zone", values = c("#A7E1BCFF", "#414083FF"), guide = guide_legend(order = 2)) +
new_scale_fill()
sintP4 = facet_plot(sintP3, panel = "depth", data=sintAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depthm), color = "grey25", size = 0.1) +
scale_fill_viridis_c("Depth (m)", option = "mako", trans = "reverse", limits = c(60, 10)) +
new_scale_fill()
sintP5 = facet_plot(sintP4, panel="K = 2", data=sintAdmixK2_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
scale_fill_manual("Lineage",values = sintKColPal[c(1,5,2,3,4,6)])
sintP6 = facet_plot( sintP5, panel="K = 3", data=sintAdmixK3_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
sintP7 = facet_plot(sintP6, panel="K = 4", data=sintAdmixK4_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
sintP8 = facet_plot(sintP7, panel="K = 5", data=sintAdmixK5_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
sintP9 = facet_plot(sintP8, panel="K = 6", data=sintAdmixK6_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
theme_tree(strip.text = element_blank(),
panel.spacing = unit(.1, "line"))
}
sintImg = image_read("../figures/icons/sint.png") %>% image_ggplot()
sintStructure = facet_widths(sintP9, widths = c(0.5, 0.025, 0.025, 0.025, 0.2, 0.1, 0.1, 0.1, 0.1))
sintStructure
Population structure
sintData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "S. intersepta") %>% dplyr::filter(!Sample %in% c("SGM007", "SGM010", "SGM013", "SGM027.2", "SGM027.3", "SGM046.1", "SGM046.3", "SGM063", "SGM152.2", "SGM152.3", "SGM169", "SGM177", "SGM179", "SGM186.1", "SGM186.3", "SGM197.2", "SGM197.3", "SGM206.2", "SGM206.3"))
sintData$depthZone = factor(sintData$depthZone)
sintData$depthZone = factor(sintData$depthZone, levels(sintData$depthZone)[c(2,1)])
sintData$bank = factor(sintData$bank)
sintData$bank = factor(sintData$bank,levels(sintData$bank)[c(5, 2, 1, 4, 3)])
sintPcadmix = read.table("../data/sint/admix/sintK2.output")
sintPcadmix %>% summarise(sum(V6),sum(V7)) ## sum(V6) sum(V7)
## 1 103.7485 99.2515sintPcadmix = sintData %>% cbind(sintPcadmix) %>% rename("cluster1" = "V6", "cluster2" = "V7") %>%dplyr::select(order(colnames(.)))
sintPcadmixClust = sintPcadmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, 0))))
sintPcadmix = sintPcadmix %>% mutate(sintPcadmixClust)
sintPcadmix$cluster = as.factor(sintPcadmix$cluster)
levels(sintPcadmix$cluster) = c("Si_Deep", "Si_1", "Admixed")
sintSiteLineages = sintPcadmix %>% dplyr::select(depthZone, cluster) %>%
group_by(depthZone) %>% count(cluster) %>% mutate(Freq = n/sum(n)) %>% apply(2, function(x) gsub("\\s+", "", x)) %>% as.data.frame()
sintSiteLineages## depthZone cluster n Freq
## 1 Shallow Si_Deep 6 0.10344828
## 2 Shallow Si_1 11 0.18965517
## 3 Shallow Admixed 41 0.70689655
## 4 Mesophotic Si_Deep 19 0.13103448
## 5 Mesophotic Si_1 12 0.08275862
## 6 Mesophotic Admixed 114 0.78620690PCA
sintCov = read.table("../data/sint/sintPcangsd.cov") %>% as.matrix()
sintPcAdmix = read.table("../data/sint/admix/sintK2.output") %>%dplyr::select(V6, V7)
sintPcAdmix %>% summarise(sum(V6),sum(V7)) ## sum(V6) sum(V7)
## 1 103.7485 99.2515sintPcAdmix = sintPcAdmix %>% rename("cluster1" = "V6", "cluster2" = "V7") %>% dplyr::select(order(colnames(.)))
sintPcaEig = eigen(sintCov)
sintPcaVar = sintPcaEig$values/sum(sintPcaEig$values)*100
head(sintPcaVar)## [1] 0.7628354 0.6680628 0.6538201 0.6487587 0.6439072 0.6427670sintPcangsd = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "S. intersepta") %>% dplyr::filter(!Sample %in% c("SGM007", "SGM010", "SGM013", "SGM027.2", "SGM027.3", "SGM046.1", "SGM046.3", "SGM063", "SGM152.2", "SGM152.3", "SGM169", "SGM177", "SGM179", "SGM186.1", "SGM186.3", "SGM197.2", "SGM197.3", "SGM206.2", "SGM206.3")) %>% dplyr::select("sample" = Sample, "bank" = bank, "depth" = depthZone, "depthm" = depthM)
sintPcangsd$sitedepth = as.factor(paste(sintPcangsd$bank, sintPcangsd$depth, sep = " "))
sintPcangsd$sitedepth = factor(sintPcangsd$sitedepth, levels(sintPcangsd$sitedepth)[c(7, 6, 3, 2, 1, 4, 5)])
sintPcangsd$bank = factor(sintPcangsd$bank)
sintPcangsd$bank = factor(sintPcangsd$bank, levels(sintPcangsd$bank)[c(5, 2, 1, 3, 4)])
sintPcangsd$depth = factor(sintPcangsd$depth)
sintPcangsd$depth = factor(sintPcangsd$depth, levels(sintPcangsd$depth)[c(2, 1)])
sintPcangsd$PC1 = sintPcaEig$vectors[,1]
sintPcangsd$PC2 = sintPcaEig$vectors[,2]
sintPcangsd$PC3 = sintPcaEig$vectors[,3]
sintPcangsd$PC4 = sintPcaEig$vectors[,4]
sintPcangsdClust = sintPcAdmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, 0))))
# pcangsdClust$clusterX = as.vector(d_clust$classification)
sintPcangsd = sintPcangsd %>% mutate(sintPcangsdClust)
sintPcangsd$cluster = as.factor(sintPcangsd$cluster)
levels(sintPcangsd$cluster) = c("Si_Deep", "Si_1", "Admixed")
sintBamsClusters = sintPcangsd %>% dplyr::select(sample, cluster) %>% dplyr::arrange(sample)
sintBamsSamples = read.delim("../data/sint/bamsNoClones", header = FALSE)
sintBamsClusters$sample = sintBamsSamples$V1
# bamsClusters = as.data.frame(bamsClusters)
write.table(x = sintBamsClusters, file = "../data/sint/sintBamsClusters", sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
sintPcangsd = merge(sintPcangsd, aggregate(cbind(mean.x = PC1, mean.y = PC2)~sitedepth, sintPcangsd, mean), by="sitedepth")sintPcangsd %>% dplyr::group_by(depth,cluster) %>% dplyr::summarize(n = n())## `summarise()` has grouped output by 'depth'. You can override using the `.groups`
## argument.## # A tibble: 6 × 3
## # Groups: depth [2]
## depth cluster n
## <fct> <fct> <int>
## 1 Shallow Si_Deep 6
## 2 Shallow Si_1 11
## 3 Shallow Admixed 41
## 4 Mesophotic Si_Deep 19
## 5 Mesophotic Si_1 12
## 6 Mesophotic Admixed 114Plot PCA
sintPcaPlot12SA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.25) +
geom_vline(xintercept = 0, color = "gray90", size = 0.25) +
geom_point(data = sintPcangsd, aes(x = PC1, y = PC2, fill = bank, shape = depth, color = bank), stroke = 0, size = 3, alpha = 0.5, show.legend = FALSE) +
geom_point(data = sintPcangsd, aes(x = mean.x, y = mean.y, fill = bank, shape = depth), color = "black", size = 3.75, alpha = 1, stroke = 0.25) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = gomPal, name = "Site") +
scale_color_manual(values = gomPal, name = "Site") +
labs(x = paste0("PC 1 (", format(round(sintPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(sintPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = guide_legend(override.aes = list(size = 4, stroke = 0.25, alpha = NA), order = 2, ncol = 1), fill = guide_legend(override.aes = list(shape = 22, size = 4, fill = gomPal[c(1,2,3,4,5)], alpha = NA), order = 1, ncol = 1)) +
theme_bw()
sintPcaPlot12S = sintPcaPlot12SA +
pcaTheme +
theme(legend.position = c(0.87, 0.76))
sintPcaPlot12S
sintPcaPlot12LA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
geom_point(data = sintPcangsd, aes(x = PC1, y = PC2, fill = cluster, shape = depth), color = "black", size = 3, alpha = 1, show.legend = TRUE) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = sintKColPal[c(1,2,7)], name = "Lineage") +
labs(x = paste0("PC 1 (", format(round(sintPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(sintPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = "none", fill = guide_legend(override.aes = list(shape = 22, size = 4, fill = sintKColPal[c(1,2,7)], alpha = NA), order = 1, ncol = 1))+
theme_bw()
sintPcaPlot12L = sintPcaPlot12LA +
pcaTheme +
theme(legend.position = c(0.9,0.86))Put all plots together
sintPcaPlots = ((sintPcaPlot12S + theme(axis.title.y = element_text(margin = ggplot2::margin(r = -20, unit = "pt")))) | sintPcaPlot12L) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
panel.background = element_rect(fill = "white"),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
sintPcaPlots
Admixture
Prepare admixture outputs for plotting
sintAdmix = sintPcangsd %>%dplyr::select(-PC1, -PC2, -PC3, -PC4, -cluster, -depthm, -mean.x, -mean.y)
sintAdmix$bank = factor(sintAdmix$bank)
sintAdmix = arrange(sintAdmix, bank, depth, -cluster1, -cluster2)
sintPopCounts = sintAdmix %>% group_by(bank, depth) %>% summarise(n = n())## `summarise()` has grouped output by 'bank'. You can override using the `.groups`
## argument.sintPopCounts## # A tibble: 7 × 3
## # Groups: bank [5]
## bank depth n
## <fct> <fct> <int>
## 1 WFGB Shallow 28
## 2 WFGB Mesophotic 30
## 3 EFGB Shallow 30
## 4 EFGB Mesophotic 28
## 5 Bright Mesophotic 30
## 6 Geyer Mesophotic 30
## 7 McGrail Mesophotic 27sintPopCountList = reshape2::melt(lapply(sintPopCounts$n,function(x){c(1:x)}))
sintAdmix$ord = sintPopCountList$value
sintAdmixMelt = melt(sintAdmix, id.vars=c("sample", "bank", "depth", "sitedepth", "ord"), variable.name="Ancestry", value.name="Fraction")
sintAdmixMelt$Ancestry = factor(sintAdmixMelt$Ancestry)
sintAdmixMelt$Ancestry = factor(sintAdmixMelt$Ancestry, levels = rev(levels(sintAdmixMelt$Ancestry)))
sintPopAnno = data.frame(x1 = c(0.5, 0.5, 0.5, 0.5, 0.5), x2 = c(30.5, 30.5, 30.5, 30.5, 30.5),
y1 = -0.1, y2 = -0.1, sample = NA, Ancestry = NA, depth = "Mesophotic",
ord = NA, Fraction = NA,
bank = c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"))
sintPopAnno$bank = factor(sintPopAnno$bank)
sintPopAnno$bank = factor(sintPopAnno$bank, levels = levels(sintPopAnno$bank)[c(5, 2, 1, 3, 4)])Make admixture plots
sintAdmixPlotA = ggplot(data = sintAdmixMelt, aes(x = ord, y = Fraction, fill = Ancestry, order = sample)) +
geom_segment(data = sintPopAnno, aes(x = x1, xend = x2, y = -.1, yend = -.1, color = bank), size = 7) +
geom_bar(stat = "identity", position = "fill", width = 1, colour = "grey25", size = 0.2) +
facet_grid(factor(depth) ~ bank, switch = "both", scales = "free_x") +
geom_text(data = (sintAdmixMelt %>% filter(depth == "Mesophotic", bank %in% c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"), sample %in% c("SGM075", "SGM127", "SGM058", "SGM208", "SGM020"), Ancestry == "cluster1")), x = 15.5, y = -.1, aes(label = bank), size = 4, color = "black") +
scale_fill_manual(values = sintKColPal) +
scale_color_manual(values = gomPal) +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
coord_cartesian(ylim = c(0.0, 1.0), clip = "off") +
theme_bw()
sintAdmixPlot = sintAdmixPlotA +
theme_bw()+
admixTheme
sintAdmixPlot
Lineage demographics
Lineage depth distribuiton
leveneTest(lm(depthm ~ cluster, data = sintPcangsd))## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.9456 0.1456
## 200sintDepthAov = welch_anova_test(depthm ~ cluster, data = subset(sintPcangsd, subset = sintPcangsd$cluster!="Admixed"))
sintDF = sintDepthAov$statistic[[1]]
sintDepthPH = games_howell_test(depthm ~ cluster, data = sintPcangsd, conf.level = 0.95) %>% as.data.frame()
sintLineageViolinA = ggplot(data = sintPcangsd, aes(x = cluster, y = depthm)) +
annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 30, ymax = Inf, fill = "black", alpha = 0.10, color = NA) +
geom_violin(aes(fill = cluster),adjust = 1, linewidth = 0, color = "black", alpha = 0.75, width = 0.9, trim = F, scale = "area") +
geom_beeswarm(aes(fill = cluster), shape = 21, size = 2, cex = 1.5, alpha = 0.5) +
geom_violin(adjust = 1, linewidth = 0.4, color = "black", alpha = 1, width = 0.9, trim = F, fill = NA, scale = "area") +
geom_boxplot(width = 0.1, color = "black", fill = "white", outlier.colour = NA, linewidth = 0.6, alpha = 0.5) +
annotate(geom = "text", x = 2.75, y =65, label = bquote(italic("F")~" = "~.(sintDF)*"; "~italic("p")~" < 0.001"), size = 3) +
scale_fill_discrete(type = sintKColPal[c(1,2,7)], name = "Lineage") +
scale_color_discrete(type = sintKColPal[c(1,2,7)], name = "Lineage") +
xlab("Lineage") +
ylab("Depth (m)") +
scale_y_reverse(breaks = seq(5, 90, 5)) +
theme_bw()
sintLineageViolin = sintLineageViolinA + violinTheme
sintLineageViolin
Symbiodiniaceae
sintData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "S. intersepta") %>% dplyr::filter(!Sample %in% c("SGM007", "SGM010", "SGM013", "SGM027.2", "SGM027.3", "SGM046.1", "SGM046.3", "SGM063", "SGM152.2", "SGM152.3", "SGM169", "SGM177", "SGM179", "SGM186.1", "SGM186.3", "SGM197.2", "SGM197.3", "SGM206.2", "SGM206.3")) %>% dplyr::select("sample" = Sample, "site" = bank, "depth" = depthZone)
sintZoox = read.delim("../data/sint/sintZooxReads", header = FALSE, check.names = FALSE)
head(sintZoox)## V1 V2
## 1 SGM001.trim.zoox.bt2.bam NA
## 2 chr1 49
## 3 chr2 75
## 4 chr3 77
## 5 chr4 107
## 6 chr5 9# Reconstruct read mapping output into dataframe usable for analysis
sintZoox$V2[is.na(sintZoox$V2)] <- as.character(sintZoox$V1[is.na(sintZoox$V2)])
sintZoox$V1 = gsub(pattern = "MGM*.", "chr", sintZoox$V1)
sintZoox$V2 = gsub(".trim.*", "", sintZoox$V2)
sintZoox = sintZoox %>% filter(sintZoox$V1 != "*")
sintZooxLst = split(sintZoox$V2, as.integer(gl(length(sintZoox$V2), 20, length(sintZoox$V2))))
sintZooxMaps = NULL
for(i in sintZooxLst){
sintZooxMaps = rbind(sintZooxMaps, data.frame(t(i)))
}
# remove tech reps
sintZooxMaps = sintZooxMaps %>% dplyr::filter(!X1 %in% c("SGM007", "SGM010", "SGM013", "SGM027.2", "SGM027.3", "SGM046.1", "SGM046.3", "SGM063", "SGM152.2", "SGM152.3", "SGM169", "SGM177", "SGM179", "SGM186.1", "SGM186.3", "SGM197.2", "SGM197.3", "SGM206.2", "SGM206.3"))
colnames(sintZooxMaps) = c("sample",sintZoox$V1[c(2:20)])
# convert characters to numeric
str(sintZooxMaps)## 'data.frame': 203 obs. of 20 variables:
## $ sample: chr "SGM001" "SGM002" "SGM003" "SGM004" ...
## $ chr1 : chr "49" "95" "217" "834" ...
## $ chr2 : chr "75" "98" "259" "1051" ...
## $ chr3 : chr "77" "107" "206" "1241" ...
## $ chr4 : chr "107" "173" "296" "1291" ...
## $ chr5 : chr "9" "19" "26" "40" ...
## $ chr6 : chr "47" "97" "266" "193" ...
## $ chr7 : chr "64" "137" "333" "239" ...
## $ chr8 : chr "76" "192" "373" "319" ...
## $ chr9 : chr "32" "114" "192" "153" ...
## $ chr10 : chr "3654" "7656" "9211" "14606" ...
## $ chr11 : chr "3936" "11339" "12614" "14468" ...
## $ chr12 : chr "4394" "12409" "14093" "16147" ...
## $ chr13 : chr "3882" "11175" "12653" "14511" ...
## $ chr14 : chr "3621" "10312" "11406" "13837" ...
## $ chr15 : chr "494" "1386" "1539" "1851" ...
## $ chr16 : chr "281" "115" "215" "936" ...
## $ chr17 : chr "93" "66" "176" "330" ...
## $ chr18 : chr "115" "61" "178" "400" ...
## $ chr19 : chr "39" "12" "33" "98" ...for(i in c(2:20)){
sintZooxMaps[,i] = as.numeric(sintZooxMaps[,i])
}
str(sintZooxMaps)## 'data.frame': 203 obs. of 20 variables:
## $ sample: chr "SGM001" "SGM002" "SGM003" "SGM004" ...
## $ chr1 : num 49 95 217 834 202 134 150 109 161 70 ...
## $ chr2 : num 75 98 259 1051 221 ...
## $ chr3 : num 77 107 206 1241 168 ...
## $ chr4 : num 107 173 296 1291 337 ...
## $ chr5 : num 9 19 26 40 29 13 18 10 8 17 ...
## $ chr6 : num 47 97 266 193 180 161 170 91 238 82 ...
## $ chr7 : num 64 137 333 239 247 180 172 163 392 104 ...
## $ chr8 : num 76 192 373 319 318 188 192 225 440 131 ...
## $ chr9 : num 32 114 192 153 131 79 84 105 202 100 ...
## $ chr10 : num 3654 7656 9211 14606 8695 ...
## $ chr11 : num 3936 11339 12614 14468 12436 ...
## $ chr12 : num 4394 12409 14093 16147 14062 ...
## $ chr13 : num 3882 11175 12653 14511 12625 ...
## $ chr14 : num 3621 10312 11406 13837 11832 ...
## $ chr15 : num 494 1386 1539 1851 1520 ...
## $ chr16 : num 281 115 215 936 64 ...
## $ chr17 : num 93 66 176 330 79 101 45 55 292 359 ...
## $ chr18 : num 115 61 178 400 52 100 39 39 234 460 ...
## $ chr19 : num 39 12 33 98 10 18 12 5 41 144 ...# collapse fake chromosomes into representative genera
sintZooxMaps$Symbiodinium = rowSums(sintZooxMaps[2:6])
sintZooxMaps$Breviolum = rowSums(sintZooxMaps[7:10])
sintZooxMaps$Cladocopium = rowSums(sintZooxMaps[11:16])
sintZooxMaps$Durusdinium = rowSums(sintZooxMaps[17:20])
# keep genera totals and turn into proportions for barplot
sintZooxMaps = sintZooxMaps[,c(1, 21:24)]
sintZooxProp = sintZooxMaps
sintZooxProp$sum = apply(sintZooxProp[, c(2:length(sintZooxProp[1,]))], 1, function(x) {
sum(x, na.rm = T)
})
sintZooxProp = cbind(sintZooxProp$sample, (sintZooxProp[, c(2:(ncol(sintZooxProp)-1))]
/ sintZooxProp$sum))
colnames(sintZooxProp)[1] = "sample"
head(sintZooxProp)## sample Symbiodinium Breviolum Cladocopium Durusdinium
## 1 SGM001 0.015062960 0.010406272 0.9494417 0.025089095
## 2 SGM002 0.008854813 0.009718698 0.9768551 0.004571387
## 3 SGM003 0.015617708 0.018106586 0.9569113 0.009364403
## 4 SGM004 0.053994791 0.010951602 0.9136834 0.021370162
## 5 SGM005 0.015140489 0.013859005 0.9677572 0.003243260
## 6 SGM006 0.012986729 0.013292814 0.9644067 0.009313715# Check that all samples total to 1
apply(sintZooxProp[, c(2:(ncol(sintZooxProp)))], 1, function(x) {
sum(x, na.rm = T)
})## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [42] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [124] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [165] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1# add sample metadata to proportions
sintSym = sintData %>% left_join(sintZooxProp)## Joining with `by = join_by(sample)`sintAdmixOrd = sintAdmix %>% dplyr::select(sample, ord)
sintPcangsdITS =sintPcangsd %>% dplyr::select(sample, cluster)
sintSym = sintSym %>% left_join(sintAdmixOrd) %>% relocate(ord,.after = depth) %>% left_join(sintPcangsdITS) %>% relocate(cluster, .after = depth)## Joining with `by = join_by(sample)`
## Joining with `by = join_by(sample)`set.seed(694)
sintSymPerm = adonis2(decostand(sintSym[, c(6:((ncol(sintSym))))], "normalize") ~ site*depth*cluster, data = sintSym, permutations = 9999, method = "bray")
sintSymPerm## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = decostand(sintSym[, c(6:((ncol(sintSym))))], "normalize") ~ site * depth * cluster, data = sintSym, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Model 20 0.22245 0.13224 1.3867 0.2549
## Residual 182 1.45973 0.86776
## Total 202 1.68217 1.00000sintSym$depth = factor(sintSym$depth)
sintSym$depth = factor(sintSym$depth, levels = levels(sintSym$depth)[c(2, 1)])
sintSym$site = factor(sintSym$site)
sintSym$site = factor(sintSym$site, levels = levels(sintSym$site)[c(5, 2, 1, 3, 4)])
#turn into melted dataframe with otustack() and remove "summ" rows
sintGssSym = otuStack(sintSym, count.columns = c(6:length(sintSym[1, ])),
condition.columns = c(1:5)) %>% filter(otu != "summ") %>% droplevels()
#check that levels are correct/ordered
levels(sintGssSym$otu)## [1] "Symbiodinium" "Breviolum" "Cladocopium" "Durusdinium"levels(sintGssSym$depth)## [1] "Shallow" "Mesophotic"levels(sintGssSym$site)## [1] "WFGB" "EFGB" "Bright" "Geyer" "McGrail"Creating Symbiodiniaceae genera relative proportion barplot
sintZooxAnno = data.frame(x1 = c(0.5, 0.5, 0.5, 0.5, 0.5), x2 = c(30.5, 30.5, 30.5, 30.5, 30.5),
y1 = -0.22, y2 = -0.22, site = c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"))
sintZooxAnno$site = factor(sintZooxAnno$site)
sintZooxAnno$site = factor(sintZooxAnno$site, levels = levels(sintZooxAnno$site)[c(5, 2, 1, 3, 4)])
sintGssSymPlot = sintGssSym %>% left_join(sintZooxAnno, by = "site") %>% left_join(sintPcangsdITS)## Joining with `by = join_by(sample, cluster)`sintGssSymPlot$ord = as.numeric(sintGssSymPlot$ord)
sintZooxPlotA = ggplot(data = sintGssSymPlot, aes(x = ord, y = count, fill = otu, order = ord)) +
geom_point(aes(x=1, y=0.5, fill = otu), shape = 22, size = 0) +
geom_bar(stat = "identity", position = "stack", colour = "grey25", width = 1, size = 0.2, show.legend = FALSE) +
xlab("Population") +
scale_fill_manual(values = colPalZoox, name = "Symbiodiniaceae \ngenus") +
geom_segment(data = (sintGssSymPlot %>% filter(depth == "Mesophotic")), aes(x = x1, xend = x2, y = y1, yend = y2, color = site), size = 7) +
scale_color_manual(values = gomPal, guide = "none") +
ggnewscale::new_scale_color() +
geom_segment(data = sintGssSymPlot, aes(x = ord-0.5, xend = ord+0.5, color = cluster), y = -.07, yend = -.07, linewidth = 4) +
scale_color_manual(values = sintKColPal[c(1,2,7)], name = "Lineage") +
coord_cartesian(ylim = c(0, 1), xlim = c(0.5, 30.5), clip = "off") +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
facet_grid(factor(depth) ~ site, drop = TRUE, scales = "free", switch = "both", space = "free") +
geom_text(data = subset(sintGssSymPlot, subset = otu == "Symbiodinium") %>% filter(sample %in% c("SGM075", "SGM127", "SGM058", "SGM208", "SGM020")), x = 15.5, y = -.205, aes(label = site), size = 4.4, color = "#000000") +
guides(fill = guide_legend(override.aes = list(size = 4), ncol = 1, order = 1), color = guide_legend(override.aes = list(size = 0.1), ncol = 1)) +
theme_bw()
sintZooxPlot = sintZooxPlotA + zooxTheme
sintZooxPlot
All S. intersepta plots
sintPlot1 = sintStructure + inset_element(sintImg, -0.1, 0.7, 0.3, 1, align_to = "full", ignore_tag = TRUE)
sintPlot2 = (sintPcaPlots + sintLineageViolin) + plot_layout(guides = "collect")
sintPlot3 = sintAdmixPlot + sintZooxPlot
sintPlots = sintPlot1 + sintPlot2 + sintPlot3 +
plot_layout(heights = c(2,1,1)) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
ggsave("../figures/figure3.png", plot = sintPlots, height = 11, width = 12, units = "in", dpi = 300)Orbicella faveolata
Dendrogram and structure
ofavBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "O. faveolata") %>% dplyr::filter(!Sample %in% c("OGM003", "OGM043", "OGM075", "OGM076", "OGM121", "OGM126", "OGM108.2", "OGM108.3", "OGM024.3", "OGM024.1", "OGM017", "OGM066.2", "OGM066.1", "OGM071", "OGM072", "OGM016", "OGM009", "OGM077", "OGM086", "OGM013", "OGM115", "OGM118", "OGM012", "OGM124", "OGM005"))# list of bams files and their populations (tech reps removed)
ofavBams$Sample <- gsub("\\.[1-3]*$", "", ofavBams$Sample)
ofavBams = ofavBams %>% dplyr::mutate(sample = Sample) %>% dplyr::select(-Sample) %>% dplyr::relocate(sample)
ofavBams$bank = factor(ofavBams$bank)
ofavBams$bank = factor(ofavBams$bank, levels = levels(mcavBams$bank))
ofavBams$depthZone = factor(ofavBams$depthZone)
ofavBams$depthZone = factor(ofavBams$depthZone, levels = levels(ofavBams$depthZone)[c(2, 1)])
ofavMa = as.matrix(read.table("../data/ofav/ofavNoClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(ofavMa) = list(ofavBams[,3],ofavBams[,3])
## admixture K = 2
#-------------------------------------
ofavK2ad = read.table("../data/ofav/admix/ofavK2.output") %>% dplyr::select(V6, V7)
ofavK2ad %>% summarise(sum(V6),sum(V7))## sum(V6) sum(V7)
## 1 67.5506 48.4494ofavAdmixK2 = ofavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(ofavK2ad) %>%
dplyr::rename("Of1" = "V6", "Of2" = "V7")
ofavAdmixK2_melt = melt(ofavAdmixK2, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 3
#-------------------------------------
ofavK3ad = read.table("../data/ofav/admix/ofavK3.output") %>% dplyr::select(V6, V7, V8)
ofavK3ad %>% summarise(sum(V6),sum(V7), sum(V8))## sum(V6) sum(V7) sum(V8)
## 1 59.9384 47.5495 8.512ofavAdmixK3 = ofavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(ofavK3ad) %>%
dplyr::rename("Of1" = "V6", "Of2" = "V7", "Of3" = "V8") %>%
dplyr::select(order(colnames(.)))
ofavAdmixK3_melt = melt(ofavAdmixK3, id = c("sample", "bank", "depth", "depthm"))
## ngsAdmix K = 4
#-------------------------------------
ofavK4ad = read.table("../data/ofav/admix/ofavK4.output") %>% dplyr::select(V6, V7, V8, V9)
ofavK4ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9))## sum(V6) sum(V7) sum(V8) sum(V9)
## 1 59.5271 28.3381 8.5588 19.5758ofavAdmixK4 = ofavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(ofavK4ad) %>%
dplyr::rename("Of1" = "V6", "Of2" = "V7", "Of3" = "V8", "Of4" = "V9") %>%
dplyr::select(order(colnames(.)))
ofavAdmixK4_melt = melt(ofavAdmixK4, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 5
#-------------------------------------
ofavK5ad = read.table("../data/ofav/admix/ofavK5.output") %>% dplyr::select(V6, V7, V8, V9, V10)
ofavK5ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10)
## 1 32.7518 37.0866 8.1358 10.7534 27.2718ofavAdmixK5 = ofavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(ofavK5ad) %>%
dplyr::rename("Of1" = "V6", "Of2" = "V7", "Of3" = "V8", "Of4" = "V9", "Of5" = "V10") #%>%
ofavAdmixK5_melt = melt(ofavAdmixK5, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 6
#-------------------------------------
ofavK6ad = read.table("../data/ofav/admix/ofavK6.output") %>% dplyr::select(V6, V7, V8, V9, V10,V11)
ofavK6ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10), sum(V11))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10) sum(V11)
## 1 32.3898 28.4808 8.0751 16.5413 22.871 7.6421ofavAdmixK6 = ofavBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(ofavK6ad) %>%
dplyr::rename("Of1" = "V6", "Of2" = "V7", "Of3" = "V8", "Of4" = "V9", "Of5" = "V10", "Of6" = "V11") #%>%
ofavAdmixK6_melt = melt(ofavAdmixK6, id = c("sample", "bank", "depth", "depthm"))
## construct figure
#-------------------------------------
# ggtree(tr) + geom_nodelab(aes(label = node), hjust = -0.5)
{
ofavTr = hclust(dist(ofavMa),"ave")
ofavP1 = ggtree(ofavTr, layout = "rectangular", size = 0.35)
ofavP2 = facet_plot(ofavP1, panel = "location", data=ofavAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = bank), color = "grey25", size = 0.1, show.legend = TRUE) +
scale_fill_manual("Site", values = gomPal, guide = guide_legend(order = 1), drop = FALSE) +
new_scale_fill()
ofavP3 = facet_plot(ofavP2, panel = "depth zone", data=ofavAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depth), color = "grey25", size = 0.1) +
scale_fill_manual("Depth zone", values = c("#A7E1BCFF", "#414083FF"), guide = guide_legend(order = 2)) +
new_scale_fill()
ofavP4 = facet_plot(ofavP3, panel = "depth", data=ofavAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depthm), color = "grey25", size = 0.1) +
scale_fill_viridis_c("Depth (m)", option = "mako", trans = "reverse", limits = c(60, 10)) +
new_scale_fill()
ofavP5 = facet_plot(ofavP4, panel="K = 2", data=ofavAdmixK2_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1) +
scale_fill_manual("Lineage", values = ofavKColPal[1:6], guide = "none")
ofavP6 = facet_plot( ofavP5, panel="K = 3", data=ofavAdmixK3_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
ofavP7 = facet_plot(ofavP6, panel="K = 4", data=ofavAdmixK4_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
ofavP8 = facet_plot(ofavP7, panel="K = 5", data=ofavAdmixK5_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
ofavP9 = facet_plot(ofavP8, panel="K = 6", data=ofavAdmixK6_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable),
stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
theme_tree(strip.text = element_blank(),
panel.spacing = unit(.1, "line"))
}
ofavImg = image_read("../figures/icons/ofav.png") %>% image_ggplot()
ofavStructure = facet_widths(ofavP9, widths = c(0.5, 0.025, 0.025, 0.025, 0.1, 0.1, 0.2, 0.1, 0.1)) #+ inset_element(ofavImg, 0.0, 0.9, 0.1, 1.0, align_to = "full", ignore_tag = TRUE)
ofavStructure
Population structure
ofavData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "O. faveolata") %>% dplyr::filter(!Sample %in% c("OGM003", "OGM043", "OGM075", "OGM076", "OGM121", "OGM126", "OGM108.2", "OGM108.3", "OGM024.3", "OGM024.1", "OGM017", "OGM066.2", "OGM066.1", "OGM071", "OGM072", "OGM016", "OGM009", "OGM077", "OGM086", "OGM013", "OGM115", "OGM118", "OGM012", "OGM124", "OGM005"))
ofavData$depthZone = factor(ofavData$depthZone)
ofavData$depthZone = factor(ofavData$depthZone, levels(ofavData$depthZone)[c(2,1)])
ofavData$bank = factor(ofavData$bank)
ofavData$bank = factor(ofavData$bank,levels(ofavData$bank)[c(5, 2, 1, 4, 3)])
ofavPcadmix = read.table("../data/ofav/admix/ofavK3.output") %>%dplyr::select(V6, V7, V8)
ofavPcadmix %>% summarise(sum(V6),sum(V7), sum(V8)) ## sum(V6) sum(V7) sum(V8)
## 1 59.9384 47.5495 8.512ofavPcadmix = ofavData %>% cbind(ofavPcadmix) %>% rename("cluster1" = "V6", "cluster2" = "V7", "cluster3" = "V8") %>%dplyr::select(order(colnames(.)))
ofavPcadmixClust = ofavPcadmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75 & cluster3 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, ifelse(cluster3 >= 0.75, 3, 0)))))
ofavPcadmix = ofavPcadmix %>% mutate(ofavPcadmixClust)
ofavPcadmix$cluster = as.factor(ofavPcadmix$cluster)
levels(ofavPcadmix$cluster) = c("Ofav1", "Ofav2", "Ofav3", "Admixed")
ofavSiteLineages = ofavPcadmix %>% dplyr::select(depthZone, cluster) %>%
group_by(depthZone) %>% count(cluster) %>% mutate(Freq = n/sum(n)) %>% apply(2, function(x) gsub("\\s+", "", x)) %>% as.data.frame()PCA
ofavCov = read.table("../data/ofav/ofavPcangsd.cov") %>% as.matrix()
ofavPcAdmix = read.table("../data/ofav/admix/ofavK3.output") %>%dplyr::select(V6, V7, V8)
ofavPcAdmix %>% summarise(sum(V6),sum(V7), sum(V8)) ## sum(V6) sum(V7) sum(V8)
## 1 59.9384 47.5495 8.512ofavPcAdmix = ofavPcAdmix %>% rename("cluster1" = "V6", "cluster2" = "V8", "cluster3" = "V7") %>% dplyr::select(order(colnames(.)))
ofavPcaEig = eigen(ofavCov)
ofavPcaVar = ofavPcaEig$values/sum(ofavPcaEig$values)*100
head(ofavPcaVar)## [1] 17.9764042 2.1801653 1.0518074 0.9772828 0.9348404 0.9206296ofavPcangsd = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "O. faveolata") %>% dplyr::filter(!Sample %in% c("OGM003", "OGM043", "OGM075", "OGM076", "OGM121", "OGM126", "OGM108.2", "OGM108.3", "OGM024.3", "OGM024.1", "OGM017", "OGM066.2", "OGM066.1", "OGM071", "OGM072", "OGM016", "OGM009", "OGM077", "OGM086", "OGM013", "OGM115", "OGM118", "OGM012", "OGM124", "OGM005")) %>% dplyr::select("sample" = Sample, "bank" = bank, "depth" = depthZone, "depthm" = depthM)
ofavPcangsd$sitedepth = as.factor(paste(ofavPcangsd$bank, ofavPcangsd$depth, sep = " "))
ofavPcangsd$sitedepth = factor(ofavPcangsd$sitedepth, levels(ofavPcangsd$sitedepth)[c(4, 3, 2, 1)])
ofavPcangsd$bank = factor(ofavPcangsd$bank)
ofavPcangsd$bank = factor(ofavPcangsd$bank, levels(ofavPcangsd$bank)[c(5, 2, 1, 4, 3)])
ofavPcangsd$depth = factor(ofavPcangsd$depth)
ofavPcangsd$depth = factor(ofavPcangsd$depth, levels(ofavPcangsd$depth)[c(2, 1)])
ofavPcangsd$PC1 = ofavPcaEig$vectors[,1]
ofavPcangsd$PC2 = ofavPcaEig$vectors[,2]
ofavPcangsd$PC3 = ofavPcaEig$vectors[,3]
ofavPcangsd$PC4 = ofavPcaEig$vectors[,4]
ofavPcangsdClust = ofavPcAdmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75 & cluster3 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, ifelse(cluster3 >= 0.75, 3, 0)))))
# pcangsdClust$clusterX = as.vector(d_clust$classification)
ofavPcangsd = ofavPcangsd %>% mutate(ofavPcangsdClust)
ofavPcangsd$cluster = as.factor(ofavPcangsd$cluster)
levels(ofavPcangsd$cluster) = c("Of_Deep1", "Of_Deep2", "Of_Shal", "Admixed")
bamsClusters = ofavPcangsd %>% dplyr::select(sample, cluster) %>% dplyr::arrange(sample)
bamsSamples = read.delim("../data/ofav/bamsNoClones", header = FALSE)
bamsClusters$sample = bamsSamples$V1
# bamsClusters = as.data.frame(bamsClusters)
write.table(x = bamsClusters, file = "../data/ofavBamsClusters", sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
ofavPcangsd = merge(ofavPcangsd, aggregate(cbind(mean.x = PC1, mean.y = PC2)~sitedepth, ofavPcangsd, mean), by="sitedepth")ofavPcangsd %>% dplyr::group_by(depth,cluster) %>% dplyr::summarize(n = n())## `summarise()` has grouped output by 'depth'. You can override using the `.groups`
## argument.## # A tibble: 5 × 3
## # Groups: depth [2]
## depth cluster n
## <fct> <fct> <int>
## 1 Shallow Of_Deep1 10
## 2 Shallow Of_Shal 47
## 3 Mesophotic Of_Deep1 51
## 4 Mesophotic Of_Deep2 7
## 5 Mesophotic Of_Shal 1Plot PCA
ofavPcaPlot12SA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.25) +
geom_vline(xintercept = 0, color = "gray90", size = 0.25) +
geom_point(data = ofavPcangsd, aes(x = PC1, y = PC2, fill = bank, shape = depth, color = bank), stroke = 0, size = 3, alpha = 0.5, show.legend = FALSE) +
geom_point(data = ofavPcangsd, aes(x = mean.x, y = mean.y, fill = bank, shape = depth), color = "black", size = 3.75, alpha = 1, stroke = 0.25) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = gomPal, name = "Site") +
scale_color_manual(values = gomPal, name = "Site") +
labs(x = paste0("PC 1 (", format(round(ofavPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(ofavPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = guide_legend(override.aes = list(size = 4, stroke = 0.25, alpha = NA), order = 2, ncol = 1), fill = guide_legend(override.aes = list(shape = 22, size = 4, fill = gomPal[c(1,2)], alpha = NA), order = 1, ncol = 1)) +
theme_bw()
ofavPcaPlot12S = ofavPcaPlot12SA +
pcaTheme +
theme(legend.position = c(0.87, 0.82))
ofavPcaPlot12S
ofavPcaPlot12LA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
geom_point(data = ofavPcangsd, aes(x = PC1, y = PC2, fill = cluster, shape = depth), color = "black", size = 3, alpha = 1, show.legend = TRUE) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = ofavKColPal[c(1,3,2)], name = "Lineage") +
labs(x = paste0("PC 1 (", format(round(ofavPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(ofavPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = "none", fill = guide_legend(override.aes = list(shape = 22, size = 4, alpha = NA, fill = ofavKColPal[c(1,3,2)]), order = 1, ncol = 1))+
theme_bw()
ofavPcaPlot12L = ofavPcaPlot12LA +
pcaTheme +
theme(legend.position = c(0.9,0.86))Put all plots together
ofavPcaPlots = ((ofavPcaPlot12S + theme(axis.title.y = element_text(margin = ggplot2::margin(r = -20, unit = "pt")))) | ofavPcaPlot12L) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
panel.background = element_rect(fill = "white"),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
ofavPcaPlots
Admixture
Prepare admixture outputs for plotting
ofavAdmix = ofavPcangsd %>%dplyr::select(-PC1, -PC2, -PC3, -PC4, -cluster, -depthm, -mean.x, -mean.y)
ofavAdmix$bank = factor(ofavAdmix$bank)
ofavAdmix = arrange(ofavAdmix, bank, depth, -cluster1, -cluster2)
ofavPopCounts = ofavAdmix %>% group_by(bank, depth) %>% summarise(n = n())## `summarise()` has grouped output by 'bank'. You can override using the `.groups`
## argument.ofavPopCounts## # A tibble: 4 × 3
## # Groups: bank [2]
## bank depth n
## <fct> <fct> <int>
## 1 WFGB Shallow 28
## 2 WFGB Mesophotic 33
## 3 EFGB Shallow 29
## 4 EFGB Mesophotic 26ofavPopCountList = reshape2::melt(lapply(ofavPopCounts$n,function(x){c(1:x)}))
ofavAdmix$ord = ofavPopCountList$value
ofavAdmixMelt = melt(ofavAdmix, id.vars=c("sample", "bank", "depth", "sitedepth", "ord"), variable.name="Ancestry", value.name="Fraction")
ofavAdmixMelt$Ancestry = factor(ofavAdmixMelt$Ancestry)
ofavAdmixMelt$Ancestry = factor(ofavAdmixMelt$Ancestry, levels = rev(levels(ofavAdmixMelt$Ancestry)))
ofavPopAnno = data.frame(x1 = c(0.5, 0.5), x2 = c(33.5, 33.5),
y1 = -0.1, y2 = -0.1, sample = NA, Ancestry = NA, depth = "Mesophotic",
ord = NA, Fraction = NA,
bank = c("WFGB", "EFGB"))
ofavPopAnno$bank = factor(ofavPopAnno$bank)
ofavPopAnno$bank = factor(ofavPopAnno$bank, levels = levels(ofavPopAnno$bank)[c(2, 1)])Make admixture plots
ofavAdmixPlotA = ggplot(data = ofavAdmixMelt, aes(x = ord, y = Fraction, fill = Ancestry, order = sample)) +
geom_segment(data = ofavPopAnno, aes(x = x1, xend = x2, y = -.1, yend = -.1, color = bank), size = 7) +
geom_bar(stat = "identity", position = "fill", width = 1, colour = "grey25", size = 0.2) +
facet_grid(factor(depth) ~ bank, switch = "both", scales = "free_x") +
geom_text(data = (ofavAdmixMelt %>% filter(depth == "Mesophotic", bank %in% c("WFGB", "EFGB"), sample %in% c(
"OGM001", "OGM073"), Ancestry == "cluster1")), x = 15.5, y = -.1, aes(label = bank), size = 4, color = "black") +
scale_fill_manual(values = ofavKColPal[c(1,3,2)]) +
scale_color_manual(values = gomPal) +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
coord_cartesian(ylim = c(0.0, 1.0), clip = "off") +
theme_bw()
ofavAdmixPlot = ofavAdmixPlotA + admixTheme
ofavAdmixPlot
Lineage demographics
Depth distribution of lineages
leveneTest(lm(depthm ~ cluster, data = ofavPcangsd))## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 8.5905 0.0003365 ***
## 113
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ofavDepthAov = welch_anova_test(depthm ~ cluster, data = subset(ofavPcangsd, subset = ofavPcangsd$cluster!="Admixed"))
ofavDepthAov## # A tibble: 1 × 7
## .y. n statistic DFn DFd p method
## * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 depthm 116 235. 2 18.8 5.15e-14 Welch ANOVAofavDF = ofavDepthAov$statistic[[1]]
ofavDepthPH = games_howell_test(depthm ~ cluster, data = subset(ofavPcangsd, subset = ofavPcangsd$cluster!="Admixed"), conf.level = 0.95) %>% as.data.frame()
ofavLineageViolinA = ggplot(data = subset(ofavPcangsd, subset = ofavPcangsd$cluster!="Admixed"), aes(x = cluster, y = depthm)) +
annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 30, ymax = Inf, fill = "black", alpha = 0.10, color = NA) +
geom_violin(aes(fill = cluster),adjust = 1, linewidth = 0, color = "black", alpha = 0.75, width = 0.9, trim = F, scale = "area") +
geom_beeswarm(aes(fill = cluster), shape = 21, size = 2, cex = 1.5, alpha = 0.5) +
geom_violin(adjust = 1, linewidth = 0.4, color = "black", alpha = 1, width = 0.9, trim = F, fill = NA, scale = "area") +
geom_boxplot(width = 0.1, color = "black", fill = "white", outlier.colour = NA, linewidth = 0.6, alpha = 0.5) +
annotate(geom = "text", x = 2.75, y =65, label = bquote(italic("F")~" = "~.(ofavDF)*"; "~italic("p")~" < 0.001"), size = 3) +
scale_fill_discrete(type = ofavKColPal[c(1,3,2)], name = "Lineage") +
scale_color_discrete(type = ofavKColPal[c(1,3,2)], name = "Lineage") +
xlab("Lineage") +
ylab("Depth (m)") +
scale_y_reverse(breaks = seq(5, 90, 5)) +
theme_bw()
ofavLineageViolin = ofavLineageViolinA + violinTheme
ofavLineageViolin
Symbiodiniaceae
ofavData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "O. faveolata") %>% dplyr::filter(!Sample %in% c("OGM003", "OGM043", "OGM075", "OGM076", "OGM121", "OGM126", "OGM108.2", "OGM108.3", "OGM024.3", "OGM024.1", "OGM017", "OGM066.2", "OGM066.1", "OGM071", "OGM072", "OGM016", "OGM009", "OGM077", "OGM086", "OGM013", "OGM115", "OGM118", "OGM012", "OGM124", "OGM005")) %>% dplyr::select("sample" = Sample, "site" = bank, "depth" = depthZone)
ofavZoox = read.delim("../data/ofav/ofavZooxReads", header = FALSE, check.names = FALSE)
head(ofavZoox)## V1 V2
## 1 OGM001.trim.zoox.bt2.bam NA
## 2 chr1 38
## 3 chr2 46
## 4 chr3 55
## 5 chr4 23
## 6 chr5 116# Reconstruct read mapping output into dataframe usable for analysis
ofavZoox$V2[is.na(ofavZoox$V2)] <- as.character(ofavZoox$V1[is.na(ofavZoox$V2)])
ofavZoox$V1 = gsub(pattern = "MGM*.", "chr", ofavZoox$V1)
ofavZoox$V2 = gsub(".trim.*", "", ofavZoox$V2)
ofavZoox = ofavZoox %>% filter(ofavZoox$V1 != "*")
ofavZooxLst = split(ofavZoox$V2, as.integer(gl(length(ofavZoox$V2), 20, length(ofavZoox$V2))))
ofavZooxMaps = NULL
for(i in ofavZooxLst){
ofavZooxMaps = rbind(ofavZooxMaps, data.frame(t(i)))
}
# remove tech reps
ofavZooxMaps = ofavZooxMaps %>% dplyr::filter(!X1 %in% c("OGM003", "OGM043", "OGM075", "OGM076", "OGM121", "OGM126", "OGM108.2", "OGM108.3", "OGM024.3", "OGM024.1", "OGM017", "OGM066.2", "OGM066.1", "OGM071", "OGM072", "OGM016", "OGM009", "OGM077", "OGM086", "OGM013", "OGM115", "OGM118", "OGM012", "OGM124", "OGM005"))
colnames(ofavZooxMaps) = c("sample",ofavZoox$V1[c(2:20)])
# convert characters to numeric
str(ofavZooxMaps)## 'data.frame': 116 obs. of 20 variables:
## $ sample: chr "OGM001" "OGM002" "OGM004" "OGM006" ...
## $ chr1 : chr "38" "30" "161" "45" ...
## $ chr2 : chr "46" "41" "354" "173" ...
## $ chr3 : chr "55" "19" "301" "214" ...
## $ chr4 : chr "23" "30" "173" "220" ...
## $ chr5 : chr "116" "88" "149" "170" ...
## $ chr6 : chr "503" "821" "28912" "5009" ...
## $ chr7 : chr "623" "1022" "32166" "5925" ...
## $ chr8 : chr "654" "1068" "31717" "6158" ...
## $ chr9 : chr "62" "102" "1448" "323" ...
## $ chr10 : chr "544" "1132" "1268" "16348" ...
## $ chr11 : chr "30" "62" "318" "515" ...
## $ chr12 : chr "54" "95" "257" "305" ...
## $ chr13 : chr "19" "61" "189" "167" ...
## $ chr14 : chr "885" "1126" "2376" "1586" ...
## $ chr15 : chr "22" "56" "86" "398" ...
## $ chr16 : chr "170" "408" "604" "3964" ...
## $ chr17 : chr "52" "113" "292" "1236" ...
## $ chr18 : chr "48" "162" "317" "1608" ...
## $ chr19 : chr "10" "45" "42" "445" ...for(i in c(2:20)){
ofavZooxMaps[,i] = as.numeric(ofavZooxMaps[,i])
}
str(ofavZooxMaps)## 'data.frame': 116 obs. of 20 variables:
## $ sample: chr "OGM001" "OGM002" "OGM004" "OGM006" ...
## $ chr1 : num 38 30 161 45 193 112 300 223 74 78 ...
## $ chr2 : num 46 41 354 173 251 319 415 365 216 96 ...
## $ chr3 : num 55 19 301 214 354 250 464 342 168 130 ...
## $ chr4 : num 23 30 173 220 284 95 191 239 129 123 ...
## $ chr5 : num 116 88 149 170 157 186 144 295 141 144 ...
## $ chr6 : num 503 821 28912 5009 25205 ...
## $ chr7 : num 623 1022 32166 5925 28709 ...
## $ chr8 : num 654 1068 31717 6158 28100 ...
## $ chr9 : num 62 102 1448 323 1185 ...
## $ chr10 : num 544 1132 1268 16348 3836 ...
## $ chr11 : num 30 62 318 515 455 294 471 292 252 195 ...
## $ chr12 : num 54 95 257 305 315 309 453 354 198 201 ...
## $ chr13 : num 19 61 189 167 186 247 378 321 138 79 ...
## $ chr14 : num 885 1126 2376 1586 2498 ...
## $ chr15 : num 22 56 86 398 181 91 181 162 106 141 ...
## $ chr16 : num 170 408 604 3964 1161 ...
## $ chr17 : num 52 113 292 1236 570 ...
## $ chr18 : num 48 162 317 1608 572 ...
## $ chr19 : num 10 45 42 445 88 64 112 118 80 97 ...# collapse fake chromosomes into representative genera
ofavZooxMaps$Symbiodinium = rowSums(ofavZooxMaps[2:6])
ofavZooxMaps$Breviolum = rowSums(ofavZooxMaps[7:10])
ofavZooxMaps$Cladocopium = rowSums(ofavZooxMaps[11:16])
ofavZooxMaps$Durusdinium = rowSums(ofavZooxMaps[17:20])
# keep genera totals and turn into proportions for barplot
ofavZooxMaps = ofavZooxMaps[,c(1, 21:24)]
ofavZooxProp = ofavZooxMaps
ofavZooxProp$sum = apply(ofavZooxProp[, c(2:length(ofavZooxProp[1,]))], 1, function(x) {
sum(x, na.rm = T)
})
ofavZooxProp = cbind(ofavZooxProp$sample, (ofavZooxProp[, c(2:(ncol(ofavZooxProp)-1))]
/ ofavZooxProp$sum))
colnames(ofavZooxProp)[1] = "sample"
head(ofavZooxProp)## sample Symbiodinium Breviolum Cladocopium Durusdinium
## 1 OGM001 0.07030855 0.4658574 0.39301973 0.07081437
## 2 OGM002 0.03209381 0.4648974 0.39068045 0.11232834
## 3 OGM004 0.01125284 0.9318995 0.04443785 0.01240977
## 4 OGM006 0.01834453 0.3886496 0.43114107 0.16186480
## 5 OGM007 0.01313892 0.8822800 0.07922587 0.02535525
## 6 OGM008 0.04146373 0.6417827 0.25279083 0.06396276# Check that all samples total to 1
apply(ofavZooxProp[, c(2:(ncol(ofavZooxProp)))], 1, function(x) {
sum(x, na.rm = T)
})## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [42] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [83] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1# add sample metadata to proportions
ofavSym = ofavData %>% left_join(ofavZooxProp)## Joining with `by = join_by(sample)`ofavAdmixOrd = ofavAdmix %>% dplyr::select(sample, ord)
ofavPcangsdITS =ofavPcangsd %>% dplyr::select(sample, cluster)
ofavSym = ofavSym %>% left_join(ofavAdmixOrd) %>% relocate(ord,.after = depth) %>% left_join(ofavPcangsdITS) %>% relocate(cluster, .after = depth)## Joining with `by = join_by(sample)`
## Joining with `by = join_by(sample)`set.seed(694)
ofavSymPerm = adonis2(decostand(ofavSym[, c(6:((ncol(ofavSym))))], "normalize") ~ site*depth*cluster, data = ofavSym, permutations = 9999, method = "bray")
ofavSymPerm## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = decostand(ofavSym[, c(6:((ncol(ofavSym))))], "normalize") ~ site * depth * cluster, data = ofavSym, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Model 8 1.2205 0.14653 2.2962 0.0155 *
## Residual 107 7.1089 0.85347
## Total 115 8.3293 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ofavSym$depth = factor(ofavSym$depth)
ofavSym$depth = factor(ofavSym$depth, levels = levels(ofavSym$depth)[c(2, 1)])
ofavSym$site = factor(ofavSym$site)
ofavSym$site = factor(ofavSym$site, levels = levels(ofavSym$site)[c(5, 2, 1, 3, 4)])
#turn into melted dataframe with otustack() and remove "summ" rows
ofavGssSym = otuStack(ofavSym, count.columns = c(6:length(ofavSym[1, ])),
condition.columns = c(1:5)) %>% filter(otu != "summ") %>% droplevels()
#check that levels are correct/ordered
levels(ofavGssSym$otu)## [1] "Symbiodinium" "Breviolum" "Cladocopium" "Durusdinium"levels(ofavGssSym$depth)## [1] "Shallow" "Mesophotic"levels(ofavGssSym$site)## [1] "WFGB" "EFGB"Symbiodiniaceae genera barplot
ofavZooxAnno = data.frame(x1 = c(0.5, 0.5), x2 = c(33.5, 33.5),
y1 = -0.22, y2 = -0.22, site = c("WFGB", "EFGB"))
ofavZooxAnno$site = factor(ofavZooxAnno$site)
ofavZooxAnno$site = factor(ofavZooxAnno$site, levels = levels(ofavZooxAnno$site)[c(5, 2, 1, 3, 4)])
ofavGssSymPlot = ofavGssSym %>% left_join(ofavZooxAnno, by = "site") %>% left_join(ofavPcangsdITS)## Joining with `by = join_by(sample, cluster)`ofavGssSymPlot$ord = as.numeric(ofavGssSymPlot$ord)
ofavZooxPlotA = ggplot(data = ofavGssSymPlot, aes(x = ord, y = count, fill = otu, order = ord)) +
geom_point(aes(x=1, y=0.5, fill = otu), shape = 22, size = 0) +
geom_bar(stat = "identity", position = "stack", colour = "grey25", width = 1, size = 0.2, show.legend = FALSE) +
xlab("Population") +
scale_fill_manual(values = colPalZoox, name = "Symbiodiniaceae \ngenus") +
geom_segment(data = (ofavGssSymPlot %>% filter(depth == "Mesophotic")), aes(x = x1, xend = x2, y = y1, yend = y2, color = site), size = 7) +
scale_color_manual(values = gomPal, guide = "none") +
ggnewscale::new_scale_color() +
geom_segment(data = ofavGssSymPlot, aes(x = ord-0.5, xend = ord+0.5, color = cluster), y = -.07, yend = -.07, linewidth = 4) +
scale_color_manual(values = ofavKColPal[c(1,3,2,7)], name = "Lineage") +
coord_cartesian(ylim = c(0, 1), xlim = c(0.5, 33.5), clip = "off") +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
facet_grid(factor(depth) ~ site, drop = TRUE, scales = "free", switch = "both", space = "free") +
geom_text(data = subset(ofavGssSymPlot, subset = otu == "Symbiodinium") %>% filter(sample %in% c("OGM001", "OGM073")), x = 15.5, y = -.205, aes(label = site), size = 4.4, color = "#000000") +
guides(fill = guide_legend(override.aes = list(size = 4), ncol = 1, order = 1), color = guide_legend(override.aes = list(size = 0.1), ncol = 1)) +
theme_bw()
ofavZooxPlot = ofavZooxPlotA + zooxTheme
ofavZooxPlot
All O. faveolata plots
ofavPlot1 = ofavStructure + inset_element(ofavImg, -0.1, 0.7, 0.3, 1, align_to = "full", ignore_tag = TRUE)
ofavPlot2 = (ofavPcaPlots + ofavLineageViolin) + plot_layout(guides = "collect")
ofavPlot3 = ofavAdmixPlot + ofavZooxPlot
ofavPlots = ofavPlot1 + ofavPlot2 + ofavPlot3 +
plot_layout(heights = c(2,1,1)) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
ggsave("../figures/figure4.png", plot = ofavPlots, height = 11, width = 12, units = "in", dpi = 300)Xestospongia muta
Dendrogram and structure
xestoBams = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "X. muta") %>% dplyr::filter(!Sample %in% c("XGM193.2", "XGM193.3", "XGM203.1", "XGM203.3", "XGM158.1", "XGM158.2", "XGM149", "XGM118", "XGM132", "XGM117", "XGM099", "XGM110", "XGM179.2", "XGM179.3", "XGM072.1", "XGM072.3", "XGM034.1", "XGM034.3", "XGM013", "XGM067", "XGM129"))# list of bams files and their populations (tech reps removed)
xestoBams$Sample <- gsub("\\.[1-3]*$", "", xestoBams$Sample)
xestoBams = xestoBams %>% dplyr::mutate(sample = Sample) %>% dplyr::select(-Sample) %>% dplyr::relocate(sample)
xestoBams$bank = factor(xestoBams$bank)
xestoBams$bank = factor(xestoBams$bank, levels = levels(xestoBams$bank)[c(5, 2, 1, 3, 4)])
xestoBams$depthZone = factor(xestoBams$depthZone)
xestoBams$depthZone = factor(xestoBams$depthZone, levels = levels(xestoBams$depthZone)[c(2, 1)])
xestoMa = as.matrix(read.table("../data/xesto/xestoNoClones.ibsMat")) # reads in IBS matrix produced by ANGSD
dimnames(xestoMa) = list(xestoBams[,3],xestoBams[,3])
## admixture K = 2
#-------------------------------------
xestoK2ad = read.table("../data/xesto/admix/xestoK2.output") %>% dplyr::select(V6, V7)
xestoK2ad %>% summarise(sum(V6),sum(V7))## sum(V6) sum(V7)
## 1 116.3335 80.6665xestoAdmixK2 = xestoBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(xestoK2ad) %>%
dplyr::rename("Xm1" = "V6", "Xm2" = "V7")
xestoAdmixK2_melt = melt(xestoAdmixK2, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 3
#-------------------------------------
xestoK3ad = read.table("../data/xesto/admix/xestoK3.output") %>% dplyr::select(V6, V7, V8)
xestoK3ad %>% summarise(sum(V6),sum(V7), sum(V8))## sum(V6) sum(V7) sum(V8)
## 1 66.6275 75.9908 54.3813xestoAdmixK3 = xestoBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(xestoK3ad) %>%
dplyr::rename("Xm1" = "V6", "Xm2" = "V7", "Xm3" = "V8") %>%
dplyr::select(order(colnames(.)))
xestoAdmixK3_melt = melt(xestoAdmixK3, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 4
#-------------------------------------
xestoK4ad = read.table("../data/xesto/admix/xestoK4.output") %>% dplyr::select(V6, V7, V8, V9)
xestoK4ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9))## sum(V6) sum(V7) sum(V8) sum(V9)
## 1 42.6422 72.7028 53.2008 28.4541xestoAdmixK4 = xestoBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(xestoK4ad) %>%
dplyr::rename("Xm1" = "V6", "Xm2" = "V9", "Xm3" = "V7", "Xm4" = "V8") %>%
dplyr::select(order(colnames(.)))
xestoAdmixK4_melt = melt(xestoAdmixK4, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 5
#-------------------------------------
xestoK5ad = read.table("../data/xesto/admix/xestoK5.output") %>% dplyr::select(V6, V7, V8, V9, V10)
xestoK5ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10)
## 1 38.9413 39.2971 51.6064 28.8596 38.2956xestoAdmixK5 = xestoBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(xestoK5ad) %>%
dplyr::rename("Xm1" = "V6", "Xm2" = "V9", "Xm3" = "V10", "Xm4" = "V7", "Xm5" = "V8") #%>%
xestoAdmixK5_melt = melt(xestoAdmixK5, id = c("sample", "bank", "depth", "depthm"))
## admixture K = 6
#-------------------------------------
xestoK6ad = read.table("../data/xesto/admix/xestoK6.output") %>% dplyr::select(V6, V7, V8, V9, V10,V11)
xestoK6ad %>% summarise(sum(V6),sum(V7), sum(V8), sum(V9), sum(V10), sum(V11))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10) sum(V11)
## 1 23.1583 33.5617 48.5514 22.3351 41.7451 27.6487xestoAdmixK6 = xestoBams %>%
dplyr::select(sample, bank, "depth" = depthZone, "depthm" = depthM) %>%
cbind(xestoK6ad) %>%
dplyr::rename("Xm1" = "V6", "Xm2" = "V9", "Xm3" = "V10", "Xm4" = "V11", "Xm5" = "V7", "Xm6" = "V8") #%>%
xestoAdmixK6_melt = melt(xestoAdmixK6, id = c("sample", "bank", "depth", "depthm"))
## construct figure
#-------------------------------------
# ggtree(tr) + geom_nodelab(aes(label = node), hjust = -0.5)
{
xestoTr = hclust(dist(xestoMa),"ave")
xestoP1 = ggtree(xestoTr, layout = "rectangular", size = 0.35)
xestoP2 = facet_plot(xestoP1, panel = "location", data=xestoAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = bank), color = "grey25", size = 0.1) +
scale_fill_manual("Site", values = gomPal, guide = guide_legend(order = 1)) +
new_scale_fill()
xestoP3 = facet_plot(xestoP2, panel = "depth zone", data=xestoAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depth), color = "grey25", size = 0.1) +
scale_fill_manual("Depth zone", values = c("#A7E1BCFF", "#414083FF"), guide = guide_legend(order = 2)) +
new_scale_fill()
xestoP4 = facet_plot(xestoP3, panel = "depth", data=xestoAdmixK2, geom = geom_tile, mapping = aes(x = 1, fill = depthm), color = "grey25", size = 0.1) +
scale_fill_viridis_c("Depth (m)", option = "mako", trans = "reverse", limits = c(60, 10)) +
new_scale_fill()
xestoP5 = facet_plot(xestoP4, panel="K = 2", data=xestoAdmixK2_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable), stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
scale_fill_manual("Lineage",values = xestoKColPal[1:6])
xestoP6 = facet_plot( xestoP5, panel="K = 3", data=xestoAdmixK3_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable), stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
xestoP7 = facet_plot(xestoP6, panel="K = 4", data=xestoAdmixK4_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable), stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
xestoP8 = facet_plot(xestoP7, panel="K = 5", data=xestoAdmixK5_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable), stat='identity', color = "grey25", size = 0.1, show.legend = FALSE)
xestoP9 = facet_plot(xestoP8, panel="K = 6", data=xestoAdmixK6_melt, geom=geom_barh, width = 1.0, mapping = aes(x = value, fill = variable), stat='identity', color = "grey25", size = 0.1, show.legend = FALSE) +
theme_tree(strip.text = element_blank(),
panel.spacing = unit(.1, "line"))
}
xestoImg = image_read("../figures/icons/xesto.png") %>% image_ggplot()
xestoStructure = facet_widths(xestoP9, widths = c(0.5, 0.025, 0.025, 0.025, 0.1, 0.1, 0.1, 0.1, 0.2))
xestoStructure
Population structure
xestoData = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "X. muta") %>% dplyr::filter(!Sample %in% c("XGM193.2", "XGM193.3", "XGM203.1", "XGM203.3", "XGM158.1", "XGM158.2", "XGM149", "XGM118", "XGM132", "XGM117", "XGM099", "XGM110", "XGM179.2", "XGM179.3", "XGM072.1", "XGM072.3", "XGM034.1", "XGM034.3", "XGM013", "XGM067", "XGM129"))
xestoData$depthZone = factor(xestoData$depthZone)
xestoData$depthZone = factor(xestoData$depthZone, levels(xestoData$depthZone)[c(2,1)])
xestoData$bank = factor(xestoData$bank)
xestoData$bank = factor(xestoData$bank,levels(xestoData$bank)[c(5, 2, 1, 4, 3)])
xestoPcadmix = read.table("../data/xesto/admix/xestoK6.output") %>%dplyr::select(V6, V7, V8, V9, V10, V11)
xestoPcadmix %>% summarise(sum(V6), sum(V7), sum(V8), sum(V9), sum(V10), sum(V11))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10) sum(V11)
## 1 23.1583 33.5617 48.5514 22.3351 41.7451 27.6487xestoPcadmix = xestoData %>% cbind(xestoPcadmix) %>% rename("cluster1" = "V6", "cluster2" = "V9", "cluster3" = "V10", "cluster4" = "V11", "cluster5" = "V7", "cluster6" = "V8") %>% dplyr::select(order(colnames(.)))
xestoPcadmixClust = xestoPcadmix %>% mutate(cluster = ifelse(cluster1 < 0.75 & cluster2 < 0.75 & cluster3 < 0.75 & cluster4 < 0.75 & cluster5 < 0.75 & cluster6 < 0.75, "NA", ifelse(cluster1 >=0.75, 1, ifelse(cluster2 >= 0.75, 2, ifelse(cluster3 >=0.75, 3, ifelse(cluster4 >= 0.75, 4, ifelse(cluster5 >=0.75, 5, ifelse(cluster6 >= 0.75, 6, 0))))))))
xestoPcadmix = xestoPcadmix %>% mutate(xestoPcadmixClust)
xestoPcadmix$cluster = as.factor(xestoPcadmix$cluster)
levels(xestoPcadmix$cluster) = c("Xm_Deep1", "Xm_Deep2", "Xm_Deep3", "Xm_Deep4", "Xm_Shal1", "Xm_Shal2", "Admixed")
xestoSiteLineages = xestoPcadmix %>% dplyr::select(depthZone, cluster) %>%
group_by(depthZone) %>% count(cluster) %>% mutate(Freq = n/sum(n)) %>% apply(2, function(x) gsub("\\s+", "", x)) %>% as.data.frame()
xestoSiteLineages## depthZone cluster n Freq
## 1 Shallow Xm_Shal1 15 0.258620690
## 2 Shallow Xm_Shal2 38 0.655172414
## 3 Shallow Admixed 5 0.086206897
## 4 Mesophotic Xm_Deep1 20 0.143884892
## 5 Mesophotic Xm_Deep2 16 0.115107914
## 6 Mesophotic Xm_Deep3 34 0.244604317
## 7 Mesophotic Xm_Deep4 16 0.115107914
## 8 Mesophotic Xm_Shal1 7 0.050359712
## 9 Mesophotic Xm_Shal2 1 0.007194245
## 10 Mesophotic Admixed 45 0.323741007PCA
xestoCov = read.table("../data/xesto/xestoPcangsd.cov") %>% as.matrix()
xestoPcAdmix = read.table("../data/xesto/admix/xestoK6.output") %>%dplyr::select(V6, V7, V8, V9, V10, V11)
xestoPcAdmix %>% summarise(sum(V6), sum(V7), sum(V8), sum(V9), sum(V10), sum(V11))## sum(V6) sum(V7) sum(V8) sum(V9) sum(V10) sum(V11)
## 1 23.1583 33.5617 48.5514 22.3351 41.7451 27.6487# xestoPcAdmix = xestoPcAdmix %>% rename("cluster1" = "V1", "cluster2" = "V2", "cluster3" = "V3", "cluster4" = "V4", "cluster5" = "V5", "cluster6" = "V6", "cluster7" = "V7") %>% dplyr::select(order(colnames(.)))
xestoPcAdmix = xestoPcAdmix %>% rename("Xm_Deep1" = "V6", "Xm_Deep2" = "V9", "Xm_Deep3" = "V10", "Xm_Deep4" = "V11", "Xm_Shal1" = "V7", "Xm_Shal2" = "V8") %>% dplyr::select(order(colnames(.)))
xestoPcaEig = eigen(xestoCov)
xestoPcaVar = xestoPcaEig$values/sum(xestoPcaEig$values)*100
head(xestoPcaVar)## [1] 7.9030025 3.5636044 2.6138877 2.1696192 1.8047309 0.8699737xestoPcangsd = read.csv("../data/nwgomMetaData.csv") %>% dplyr::filter(species == "X. muta") %>% dplyr::filter(!Sample %in% c("XGM193.2", "XGM193.3", "XGM203.1", "XGM203.3", "XGM158.1", "XGM158.2", "XGM149", "XGM118", "XGM132", "XGM117", "XGM099", "XGM110", "XGM179.2", "XGM179.3", "XGM072.1", "XGM072.3", "XGM034.1", "XGM034.3", "XGM013", "XGM067", "XGM129")) %>% dplyr::select("sample" = Sample, "bank" = bank, "depth" = depthZone, "depthm" = depthM)
xestoPcangsd$sitedepth = as.factor(paste(xestoPcangsd$bank, xestoPcangsd$depth, sep = " "))
xestoPcangsd$sitedepth = factor(xestoPcangsd$sitedepth, levels(xestoPcangsd$sitedepth)[c(7, 6, 3, 2, 1, 4, 5)])
xestoPcangsd$bank = factor(xestoPcangsd$bank)
xestoPcangsd$bank = factor(xestoPcangsd$bank, levels(xestoPcangsd$bank)[c(5, 2, 1, 3, 4)])
xestoPcangsd$depth = factor(xestoPcangsd$depth)
xestoPcangsd$depth = factor(xestoPcangsd$depth, levels(xestoPcangsd$depth)[c(2, 1)])
xestoPcangsd$PC1 = xestoPcaEig$vectors[,1]
xestoPcangsd$PC2 = xestoPcaEig$vectors[,2]
xestoPcangsd$PC3 = xestoPcaEig$vectors[,3]
xestoPcangsd$PC4 = xestoPcaEig$vectors[,4]
xestoPcangsdClust = xestoPcAdmix %>% mutate(cluster = ifelse(Xm_Deep1 < 0.75 & Xm_Deep2 < 0.75 & Xm_Deep3 < 0.75 & Xm_Deep4 < 0.75 & Xm_Shal1 < 0.75 & Xm_Shal2 < 0.75, "Admixed", ifelse(Xm_Deep1 >=0.75, "Xm_Deep1", ifelse(Xm_Deep2 >= 0.75, "Xm_Deep2", ifelse(Xm_Deep3 >=0.75, "Xm_Deep3", ifelse(Xm_Deep4 >= 0.75, "Xm_Deep4", ifelse(Xm_Shal1 >=0.75, "Xm_Shal1", ifelse(Xm_Shal2 >= 0.75, "Xm_Shal2", 0))))))))
# pcangsdClust$clusterX = as.vector(d_clust$classification)
xestoPcangsd = xestoPcangsd %>% mutate(xestoPcangsdClust)
xestoPcangsd$cluster = as.factor(xestoPcangsd$cluster)
xestoPcangsd$cluster = factor(xestoPcangsd$cluster, levels = levels(xestoPcangsd$cluster)[c(2:7,1)])
xestoBamsClusters = xestoPcangsd %>% dplyr::select(sample, cluster) %>% dplyr::arrange(sample)
xestoBamsSamples = read.delim("../data/xesto/bamsNoClones", header = FALSE)
xestoBamsClusters$sample = xestoBamsSamples$V1
# bamsClusters = as.data.frame(bamsClusters)
write.table(x = xestoBamsClusters, file = "../data/xesto/xestoBamsClusters", sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
xestoPcangsd = merge(xestoPcangsd, aggregate(cbind(mean.x = PC1, mean.y = PC2)~sitedepth, xestoPcangsd, mean), by="sitedepth")xestoPcangsd %>% dplyr::group_by(depth,cluster) %>% dplyr::summarize(n = n())## `summarise()` has grouped output by 'depth'. You can override using the `.groups`
## argument.## # A tibble: 10 × 3
## # Groups: depth [2]
## depth cluster n
## <fct> <fct> <int>
## 1 Shallow Xm_Shal1 15
## 2 Shallow Xm_Shal2 38
## 3 Shallow Admixed 5
## 4 Mesophotic Xm_Deep1 20
## 5 Mesophotic Xm_Deep2 16
## 6 Mesophotic Xm_Deep3 34
## 7 Mesophotic Xm_Deep4 16
## 8 Mesophotic Xm_Shal1 7
## 9 Mesophotic Xm_Shal2 1
## 10 Mesophotic Admixed 45Plot PCA
xestoPcaPlot12SA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.25) +
geom_vline(xintercept = 0, color = "gray90", size = 0.25) +
geom_point(data = xestoPcangsd, aes(x = PC1, y = PC2, fill = bank, shape = depth, color = bank), stroke = 0, size = 3, alpha = 0.5, show.legend = FALSE) +
geom_point(data = xestoPcangsd, aes(x = mean.x, y = mean.y, fill = bank, shape = depth), color = "black", size = 3.75, alpha = 1, stroke = 0.25) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = gomPal, name = "Site") +
scale_color_manual(values = gomPal, name = "Site") +
labs(x = paste0("PC 1 (", format(round(xestoPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(xestoPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = guide_legend(override.aes = list(size = 4, stroke = 0.25, alpha = NA), order = 2, ncol = 1), fill = guide_legend(override.aes = list(shape = 22, size = 4, fill = gomPal[c(1,2,3,4,5)], alpha = NA), order = 1, ncol = 1)) +
theme_bw()
xestoPcaPlot12S = xestoPcaPlot12SA +
pcaTheme +
theme(legend.position = c(0.1, 0.24))
xestoPcaPlot12S
xestoPcaPlot12LA = ggplot() +
geom_hline(yintercept = 0, color = "gray90", size = 0.5) +
geom_vline(xintercept = 0, color = "gray90", size = 0.5) +
geom_point(data = xestoPcangsd, aes(x = PC1, y = PC2, fill = cluster, shape = depth), color = "black", size = 3, alpha = 1, show.legend = TRUE) +
scale_shape_manual(values = c(21, 23), name = "Depth Zone") +
scale_fill_manual(values = xestoKColPal, name = "Lineage") +
labs(x = paste0("PC 1 (", format(round(xestoPcaVar[1], 1), nsmall = 1)," %)"), y = paste0("PC 2 (", format(round(xestoPcaVar[2], 1), nsmall = 1), " %)")) +
guides(shape = "none", fill = guide_legend(override.aes = list(shape = 22, size = 4, alpha = NA), order = 1, ncol = 1))+
theme_bw()
xestoPcaPlot12L = xestoPcaPlot12LA +
pcaTheme +
theme(legend.position = c(0.9, 0.2))Put all plots together
xestoPcaPlots = ((xestoPcaPlot12S + theme(axis.title.y = element_text(margin = ggplot2::margin(r = -20, unit = "pt")))) | xestoPcaPlot12L) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
panel.background = element_rect(fill = "white"),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
xestoPcaPlots
Admixture
Prepare admixture outputs for plotting
xestoAdmix = xestoPcangsd %>%dplyr::select(-PC1, -PC2, -PC3, -PC4, -cluster, -depthm, -mean.x, -mean.y)
xestoAdmix$bank = factor(xestoAdmix$bank)
xestoAdmix = arrange(xestoAdmix, bank, depth, -Xm_Deep1, -Xm_Deep4, -Xm_Shal2)
xestoPopCounts = xestoAdmix %>% group_by(bank, depth) %>% summarise(n = n())## `summarise()` has grouped output by 'bank'. You can override using the `.groups`
## argument.xestoPopCounts## # A tibble: 7 × 3
## # Groups: bank [5]
## bank depth n
## <fct> <fct> <int>
## 1 WFGB Shallow 29
## 2 WFGB Mesophotic 29
## 3 EFGB Shallow 29
## 4 EFGB Mesophotic 26
## 5 Bright Mesophotic 25
## 6 Geyer Mesophotic 30
## 7 McGrail Mesophotic 29xestoPopCountList = reshape2::melt(lapply(xestoPopCounts$n,function(x){c(1:x)}))
xestoAdmix$ord = xestoPopCountList$value
xestoAdmixMelt = melt(xestoAdmix, id.vars=c("sample", "bank", "depth", "sitedepth", "ord"), variable.name="Ancestry", value.name="Fraction")
xestoAdmixMelt$Ancestry = factor(xestoAdmixMelt$Ancestry)
xestoAdmixMelt$Ancestry = factor(xestoAdmixMelt$Ancestry, levels = rev(levels(xestoAdmixMelt$Ancestry)))
xestoPopAnno = data.frame(x1 = c(0.5, 0.5, 0.5, 0.5, 0.5), x2 = c(30.5, 30.5, 30.5, 30.5, 30.5),
y1 = -0.1, y2 = -0.1, sample = NA, Ancestry = NA, depth = "Mesophotic",
ord = NA, Fraction = NA,
bank = c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"))
xestoPopAnno$bank = factor(xestoPopAnno$bank)
xestoPopAnno$bank = factor(xestoPopAnno$bank, levels = levels(xestoPopAnno$bank)[c(5, 2, 1, 3, 4)])Make admixture plots
xestoAdmixPlotA = ggplot(data = xestoAdmixMelt, aes(x = ord, y = Fraction, fill = Ancestry, order = sample)) +
geom_segment(data = xestoPopAnno, aes(x = x1, xend = x2, y = -.1, yend = -.1, color = bank), size = 7) +
geom_bar(stat = "identity", position = "fill", width = 1, colour = "grey25", size = 0.2) +
facet_grid(factor(depth) ~ bank, switch = "both", scales = "free_x") +
geom_text(data = (xestoAdmixMelt %>% filter(depth == "Mesophotic", bank %in% c("WFGB", "EFGB", "Bright", "Geyer", "McGrail"), sample %in% c("XGM101", "XGM135", "XGM078", "XGM177", "XGM058"), Ancestry == "Xm_Deep1")), x = 15.5, y = -.1, aes(label = bank), size = 4, color = "black") +
scale_fill_manual(values = xestoKColPal[c(1:6)]) +
scale_color_manual(values = gomPal) +
scale_x_discrete(expand = c(0.005, 0.005)) +
scale_y_continuous(expand = c(0.001, 0.001)) +
coord_cartesian(ylim = c(0.0, 1.0), clip = "off") +
theme_bw()
xestoAdmixPlot = xestoAdmixPlotA +
theme_bw()+
admixTheme
xestoAdmixPlot
Lineage demographics
Lineage depth distribution
leveneTest(lm(depthm ~ cluster, data = xestoPcangsd))## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 6 5.0923 0.00007161 ***
## 190
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1xestoDepthAov = welch_anova_test(depthm ~ cluster, data = subset(xestoPcangsd, subset = xestoPcangsd$cluster!="Admixed"))
xestoDF = xestoDepthAov$statistic[[1]]
xestoDepthPH = games_howell_test(depthm ~ cluster, data = xestoPcangsd, conf.level = 0.95) %>% as.data.frame()
xestoLineageViolinA = ggplot(data = xestoPcangsd, aes(x = cluster, y = depthm)) +
annotate(geom = "rect", xmin = -Inf, xmax = Inf, ymin = 30, ymax = Inf, fill = "black", alpha = 0.10, color = NA) +
geom_violin(aes(fill = cluster),adjust = 1, linewidth = 0, color = "black", alpha = 0.75, width = 0.7, trim = F, scale = "width") +
geom_beeswarm(aes(fill = cluster), shape = 21, size = 2, cex = 1.5, alpha = 0.5) +
geom_violin(adjust = 1, linewidth = 0.4, color = "black", alpha = 1, width = 0.7, trim = F, fill = NA, scale = "width") +
geom_boxplot(width = 0.2, color = "black", fill = "white", outlier.colour = NA, linewidth = 0.6, alpha = 0.5) +
annotate(geom = "text", x = 2.75, y =65, label = bquote(italic("F")~" = "~.(xestoDF)*"; "~italic("p")~" < 0.001"), size = 3) +
scale_fill_discrete(type = xestoKColPal, name = "Lineage") +
scale_color_discrete(type = xestoKColPal, name = "Lineage") +
xlab("Lineage") +
ylab("Depth (m)") +
scale_y_reverse(breaks = seq(5, 90, 5)) +
theme_bw()
xestoLineageViolin = xestoLineageViolinA + violinTheme + theme(axis.text.x = element_text(color = "black", size = 10, angle = 45, hjust = 0.9))
xestoLineageViolin
All X. muta plots
xestoPlot1 = xestoStructure + inset_element(xestoImg, -0.1, 0.7, 0.3, 1, align_to = "full", ignore_tag = TRUE)
xestoPlot2 = xestoPcaPlots + plot_layout(guides = "collect")
xestoPlot3 = xestoAdmixPlot + xestoLineageViolin
xestoPlots = xestoPlot1 + xestoPlot2 + xestoPlot3 +
plot_layout(heights = c(2,1,1)) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(size = 18),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
ggsave("../figures/figure5.png", plot = xestoPlots, height = 10, width = 11, units = "in", dpi = 300)Other plots
gomTrees = mcavStructure + inset_element(mcavImg, 0.02, 0.85, 0.17, 1.0, align_to = "full", ignore_tag = TRUE) +
ofavStructure + inset_element(ofavImg, 0.02, 0.85, 0.17, 1.0, align_to = "full", ignore_tag = TRUE) +
sintStructure + inset_element(sintImg, 0.02, 0.85, 0.17, 1.0, align_to = "full", ignore_tag = TRUE) +
xestoStructure + inset_element(xestoImg, 0.02, 0.85, 0.17, 1.0, align_to = "full", ignore_tag = TRUE) +
plot_annotation(tag_levels = c("A")) +
theme(plot.tag = element_text(size = 18),
legend.spacing = unit(-5, "pt"),
legend.key = element_blank(),
legend.background = element_blank())
# gomTrees
ggsave("../figures/extras/nwgom_trees.png", plot = gomTrees, dpi = 300, height = 10, width = 12, units = "in")Saving R data
save.image("nwgom.RData")